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ABSTRACT 


^Spline  functions,  because  of  their  highly  desirable  interpolating  and  approximat¬ 
ing  characteristics,  are  used  as  a  potential  alternative  to  the  conventional  pulse 
approximation  method  in  digital  image  processing.  For  uniformly  spaced  knots, 
a  class  of  spline  functions  called  B- splines  has  the  useful  properties  of  shift 
invariance,  positiveness,  and  convolutional  and  local  basis  properties.  These 
properties  are  exploited  in  image  processing  for  linear  incoherent  imaging  systems.”*' 

The  problem  of  image  degradation  in  a  linear  imaging  system  is  described  by  a 
superposition  integral.  For  simulation  of  degradation  and  restoration  by  means  of 
a  digital  computer,  the  continuous  imaging  model  must  be  discretized.  Thus,  a 
theoretical  and  experimented  study  of  quadrature  formulae,  particularly  monospline 
and  best  quadrature  formulae  in  the  sense  of  Sard,  is  presented.  It  is  shown  that 
a  good  choice  of  degree  for  a  monospline  highly  depends  on  the  frequency  content 
of  the  integrand,  and  in  most  cases,  a  cubic  monospline  generates  less  error  than 
the  pulse  approximation  method  and  Newton- Cotes  quadrature  formulae. 

In  space-invariant  imaging  systems,  the  object  and  point- spread  function  are 
represented  by  B- splines  of  degrees  m  and  n.  Exploiting  the  convolutional 
property,  the  deterministic  part  of  the  blurred  image  is  a  spline  function  of  degree 
m  +  n  +  1.  A  minimum  norm  principle  leading  to  pseudo-inversion  is  used  for  the 
restoration  of  space-variant  degradations  and  underdetermined  and  overdetermined 
models.  Space- variant  point-spread  functions  that  describe  astigmatism  and 
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curvature-of-field  are  derived  and  coordinate 
transformations  are  applied  to  reduce  the  dimen¬ 
sionality.  The  singular-value-decomposition 
technique  is  used  for  solution  of  the  simplified 
equations. 

For  noisy  blurred  images,  a  controllable 
smoothing  criteria  based  on  the  locally  variable 
statistics  and  minimization  of  the  second  deriva¬ 
tive  is  defined,  and  the  corresponding  filter, 
applicable  to  both  space- variant  and  space - 
invariant  degradations,  is  obtained.  The  para¬ 
meters  of  the  filter  determine  the  local  smooth¬ 
ing  window  and  overall  extent  of  smoothing,  and 
thus  the  trade-off  between  resolution  and  smooth¬ 
ing  is  controllable  in  a  spatially  nonstationary 
manner.  Since  the  matrices  of  this  filter  are 
banded  circulant  or  Toeplitz,  efficient  algorithms 
are  used  for  matrix  manipulations. 


14.  Key  Words:  Digital  Image  Processing, 

Monospline  Quadrature  Formulae, 
Image  Restoration,  Spline  Func¬ 
tions,  Smoothing  Splines. 
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ABSTRAC  T 


Spline  functions,  because  of  their  highly  desirable  interpolating 
and  approximating  characteristics,  are  used  as  a  potential  alterna¬ 
tive  to  the  conventional  pulse  approximation  method  in  digital  image 
processing.  For  uniformly  spaced  knots,  a  class  of  spline  functions 
called  B-splines  has  the  useful  properties  of  shift  invariance,  posi¬ 
tiveness,  and  convolutional  and  local  basis  properties.  These 
properties  are  exploited  in  image  processing  for  linear  incoherent 
imaging  systems. 

The  problem  of  image  degradation  in  a  linear  imaging  system  is 
described  by  a  superposition  integral.  For  simulation  of  degrada¬ 
tion  and  restoration  by  means  of  a  digital  computer,  the  continuous 
imaging  model  must  be  discretized.  Thus,  a  theoretical  and  experi¬ 
mental  study  of  quadrature  formulae,  particularly  monospline  and 
beat  quadrature  formulae  in  the  sense  of  Sard,  is  presented.  It  is 
shown  that  a  good  choice  of  degree  for  a  monospline  highly  depends 
on  the  frequency  content  of  the  integrand,  and  in  most  cases,  a 
cubic  monospline  generates  less  error  than  the  pulse  approximation 
method  and  Newton-Cotes  quadrature  formulae. 

In  space-invariant  imaging  systems,  the  object  and  point-spread 
function  are  represented  by  B-splines  of  degrees  m  and  n.  Kxploi- 
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ting  the  convolutional  property,  the  deterministic  part  of  the  blurred 
image  is  a  spline  function  of  degree  m+n  +  1.  A  minimum  norm  prin¬ 
ciple  leading  to  pseudo-inversion  is  used  for  the  restoration  of  space- 
variant  degradations  and  underdetermined  and  overdetermined  mod¬ 
els.  Space -variant  point-spread  functions  that  describe  astigmat¬ 
ism  and  curvature-of -field  are  derived  and  coordinate  transforma¬ 
tions  are  applied  to  reduce  the  dimensionality.  The  singular -value  - 
decomposition  technique  is  used  for  solution  of  the  simplified  equa¬ 
tions. 

For  noisy  blurred  images,  a  controllable  smoothing  criteria 
based  on  the  locally  variable  statistics  and  minimization  of  the 
second  derivative  is  defined,  and  the  corresponding  filter,  appli¬ 
cable  to  both  space-variant  and  space-invariant  degradations,  is 
obtained.  The  parameters  of  the  filter  determine  the  local  smooth¬ 
ing  window  and  overall  extent  of  smoothing,  and  thus  the  trade-off 
between  resolution  and  smoothing  is  controllable  in  a  spatially  non¬ 
stationary  manner.  Since  the  matrices  of  this  filter  are  banded  cir- 
culant  or  Toeplitz,  efficient  algorithms  are  used  for  matrix  manipu¬ 
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Chapter  1 


IN  I'RODUC  HON 

The  objective  of  image  restoration  is  the  reconstruction  of  a 
recorded  image  towards  an  ideal  object  by  inversion  of  the  degrading 
phenomena.  These  phenomena  include  such  imperfect  imaging  cir¬ 
cumstances  as  defocus,  motion  blur,  optical  aberrations,  and  noise 
[n  ,  [2".  The  pioneers  of  this  field  in  the  modern  sense  were 
Marechal  and  his  co-workers  [3"|  who  recognized  in  the  1950’s  the 
potential  of  optical  spatial  filtering  for  restoring  blurred  photographs. 
Their  success  stimulated  others  to  study  image  restoration  by 
opti  cal  compensation  of  the  degradations.  However,  it  was  the 
versatility  of  digital  computers  and  the  space  program  of  the  sixties 
with  its  need  for  high  quality  imagery  that  provided  the  necessary 
means  and  motivation  for  the  development  of  the  field.  With  digital 
processing  it  is  possible  to  overcome  many  limitations  of  optical 
filtering  and  to  explore  new  approaches  which  have  no  conceivable 
optical  counterparts. 

Restoration  techniques  require  some  knowledge  concerning  the 
degradation  phenomena,  and  this  knowledge  may  come  from  an 
analytical  model,  statistical  model,  or  other  a  priori  information  of 
the  imaging  system.  Thus  considerable  emphasis  must  he  placed  on 
the  sources  and  models  of  degradation.  In  general,  an  exact 


degradation  model  is  too  complicated  to  be  used.  However,  for  many 


cases  of  practical  interest,  a  quite  accurate  model  is  given  by  a 
linear  smoothing  operation  due  to  the  optical  imperfection  followed 
by  the  addition  of  noise  [4]  . 

The  earlier  restoration  techniques,  mostly  optically  oriented, 
attempted  to  remove  the  degradation  by  inverse  filtering  [5*1  .  Using 
the  Fourier  transform  properties  of  lenses,  the  Fourier  transform 
of  the  degraded  image  is  simply  multiplied  by  the  inverse  of  the 
Fourier  transform  of  the  blurring  function.  This  method  is  not  without 
limitations  and  shortcomings.  First,  in  many  practical  cases  such  as 
motion  blur  and  defocusing,  the  Fourier  transform  of  the  blur  function 
has  zeros  at  spatial  frequencies  within  the  range  of  interest,  and 
since  the  inverse  of  zero  is  undefined,  the  method  breaks  down  for 
such  cases.  Second,  for  noisy  images,  this  method  enhances  the 
high  frequency  component  of  the  noise.  Various  modifications  have 
been  suggested  to  overcome  these  drawbacks,  but  all  of  them  are 
ad  hoc  and  intuitive  [3l  ,  [5],  [h]  ,  [7]  ,  [8]  .  In  spite  of  all  the  limita¬ 
tions,  inverse  filtering  can  yield  reasonably  good  results  where  noise 
is  not  the  limiting  degradation. 

The  minimum  mean-squared  error  (MSE)  criterion  has  been 
used  as  an  objective  criteria  for  restoration  of  noisy  images. 

Assuming  the  object  and  noise  to  be  uncorrelated  random  processes 
with  a  known  blur  function,  Helstrom  [9]  has  proposed  a  filter  for 
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image  restoration  based  on  minimum  MSE  principle.  Phis  filter  is 
the  same  as  the  classical  Wiener  filter  which  was  developed  in  the 
1940’s  in  the  field  of  signal  processing.  For  an  unknown  blur 
function,  Slepian  [lo]  has  solved  the  same  problem  assuming  the 
blur  function  to  be  a  random  process.  Utilizing  the  transform 
properties  of  imaging  systems,  Pratt  [l  1  ]  has  introduced  generalized 
Wiener  filtering  with  improved  computational  efficiency.  Habibi  [12] 
has  shown  that  a  lower  triangular  transformation  yields  an  efficient 
suboptimal  Wiener  filter.  The  Fourier  transform  properties  of  the 
circulant  matrices  has  been  used  to  develop  a  computationally  fast 
algorithm  for  solving  the  Wiener  filter  [l3]  . 

The  Wiener  filter  has  limitations  and  shortcomings.  The 
minimum  MSE  principle,  which  is  the  objective  criteria  of  a  Wiener 
filter,  is  suspect  in  image  restoration.  It  is  well  known  that  the 
human  visual  system  demands  a  more  accurate  reproduction  of 
regions  where  the  intensity  changes  rapidly  than  of  the  regions  with 
little  change,  and  the  sensitivity  of  the  eye  to  a  given  error  in  intensity 
depends  strongly  upon  the  intensity  [4l  .  The  minimum  MSEi  weights 
errors  equally  regardless  of  the  intensity  and  its  gradient.  More¬ 
over,  a  Wiener  filter  requires  extensive  a  priori  information, 
namely,  the  blur  function  and  detailed  knowledge  of  object  and  noise 
covariance  functions.  Finally,  since  the  Wiener  filter  is  derived  by 


the  Fourier  transform  properties  of  space -invariant  degradations  and 


stationary  assumptions  of  the  object  and  noise,  the  filter  is  not 


applicable  to  space-variant  degrad'  ions  and  non  -  stati  onary  objects 

Constrained  restoration  has  been  introduced  as  an  alternative  to 
overcome  some  short-comings  of  the  Wiener  filter.  Hunt  [l 4 "1  has 
proposed  a  constrained  least  square  filter,  in  which  by  judicious 
choice  of  some  variables  one  can  minimize  higher  order  derivatives, 
eye  model  effects,  or  even  achieve  the  Wiener  filter.  Stockham  and 
Cole  [l5*|  have  suggested  a  geometrical  mean  filter  between  the 
inverse  filter  and  Wiener  filter.  Utilizing  linear  equality  and  in¬ 
equality  constraints  has  led  to  constrained  restoration  techniques  [16"I. 
The  non-negative  nature  of  image  intensity  has  been  the  leading  factor 
in  some  restoration  techniques  [l7*J.  For  unknown  blur  functions, 
the  concept  of  homomorphic  systems  [l  fj  1  has  been  employed  to 
estimate  the  point-spread  function  from  the  degraded  image  by  taking 
averages  of  image  segments  in  the  log-spectral  domain  [19  i  .  A 
detailed  comparison  of  these  restoration  techniques  is  given  by  Hunt 

fco'’ . 

For  space-variant  degradations,  the  problem  of  image  restoration 
is  much  more  difficult  because  Fourier  techniques  cannot  be  used. 
Generally,  there  are  twice  as  many  independent  variables  in  a 
space-variant  system  as  in  a  space-invariant  system,  and  this 
increased  dimensionality  is  the  major  analytical  and  computational 
difficulty.  A  method,  analogous  to  Fourier  techniques,  has  been 
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pri-sfiv. ifl  in  Lcrms  of  the  degrading  system  eigenfunction  [2)  , 

however,  it  i.s  not  known  how  to  find  a  complete  set  of  eigenfunctions 
or  even  if  they  exist.  Sawchuk  [22]  has  sho^'n  that  for  certain  space- 
variant  systems  the  degradation  can  be  transformed  to  be  space- 
invariant  by  an  appropriate  selection  of  coordinates, 

1’he  following  is  an  outline  of  this  dissertation  and  a  summary 
of  the  contributions. 

In  Chapter  2  background  on  the  problem  of  image  degradation  and 
restoration  in  a  continuous  model  is  discussed.  The  mathematical 
representation  of  this  model,  inverse  filtering  and  the  Wiener  filter 
are  studied  briefly. 

Chapter  3  is  devoted  to  discrete  representations  of  the  continuous 
model  for  implementation  on  a  digital  computer.  The  pulse  approxi¬ 
mation  method  has  been  the  simplest  and  the  most  common  method  in 
digital  image  processing,  however,  the  accuracy  of  this  technique  is 
suspect  in  image  discretization.  It  is  shown  that  numerical  analysis 
techniques,  particularly  monospline  quadrature  formulae,  lead  to  a 
more  accurate  discrete  model.  The  results  are  compared  with 
extreme  cases,  namely,  the  pulse  approximation  method  and  the 
Newton-Cotes  quadrature  formulae.  B -splines,  because  of  their 
desirable  characteristics  and  the  useful  properties  of  shift  invariance, 
positiveness,  and  their  convolutional  and  local  basis  properties,  are 
studied  and  suggested  for  discrete  representation  of  the  continuous 
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model. 


The  restoration  of  noiseless  images  are  presented  in  Chapter  4. 
For  space-invariant  imaging  systems,  the  object  and  point-spread 
function  are  represented  by  B-splines  of  degrees  m  and  n.  I'he 
degree  of  B-spline  must  be  selected  with  respect  to  the  continuity 
and  frequency  content  of  the  approximated  function.  Exploiting  the 
convolutional  property,  the  blurred  image  is  a  B-spline  of  degree 
m  +  nll.  It  is  shown  that  B-splinc  produces  a  better  quality  restora¬ 
tion  than  the  conventional  pulse  approximation  method.  Pseudo¬ 
inversion  based  on  the  minimum  norm  principle  is  used  for  the 
restoration  of  space-variant  degradations,  overdetermined  models 
and  underdetermined  models.  With  a  linear  incoherent  system,  the 
space-variant  point-spread  functions  that  describe  imaging  in  the 
presence  of  astigmatism  and  cur''ature-of-field  are  derived  and 
coordinate  transformations  are  applied  to  reduce  the  dimensionality. 
The  singular-value-decomposition  techniques  analogous  to  inverse 
Fourier  filtering  are  used  for  pseudo-inverse  solution  of  the  simplified 
equations. 

Image  restoration  by  spline  functions  in  the  presence  of  noise  is 
covered  in  Chapter  5,  A  controllable  smoothing  criteria  based  on  the 
locally  variable  statistics  and  minimization  of  the  second  derivative  is 
defined,  and  the  corresponding  filter,  applicable  to  both  space- 
invariant  and  space -variant  degradations,  is  obtained.  The 
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parameters  of  the  filter  determine  the  local  smoothing  window  and 
overall  extent  of  smoothing,  and  thus  the  tradeoff  between  resolution 
and  smoothing  is  controllable  in  a  spatially  non  -  stationary  manner. 

The  interesting  properties  of  this  filter  has  made  it  capable  of 
restoring  signal-dependent  noisy  images,  and  it  has  been  successfully 
applied  for  filtering  images  degraded  by  film-grain  noise.  Since  the 
matrices  of  this  filter  are  banded,  circulant  or  Toeplitz,  efficient 
algorithms  are  used  for  matrix  manipulations. 

Finally,  conclusions  and  recommendations  for  further  research 
are  given  in  Chapter  6. 
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Chapter  2 


I  MACK  RESTORA  HON  IN  A  CONTINUOUS  MODEL 

In  this  chapter,  the  problem  of  image  degradation  and  restoration 
in  a  continuous  model  is  discussed.  Section  2.  1  presents  a  mathe¬ 
matical  model  for  a  linear  imaging  system.  In  sections  2.2  and  2.1 
respectively,  restoration  techniques  for  noiseless  and  noisy  images 
are  discussed. 

2.  1  Degradation  in  a  Linear  Imaging  System 

Lot  g(x,  y)  be  the  image  of  an  object  f(r,h)  which  has  been 
degraded  by  the  linear  operator  h(x,  y;",  Ti)  such  that 

» 

g(x,y)  =  Jj  htx,  y;c,  "P )ff  ”,  Tl)de  dH  +  n(x,  y)  .  (2.1) 

-  CO 

The  first  source  of  degradation,  represented  by  h(*  ),  is  known  as  the 
impulse  response  or  point  spread  function  (PSF)  of  the  imaging  system. 
Physically,  h(.  )  is  assumed  to  be  the  image  of  a  point  source  of  light 
located  at  ( ^ ,  T)  >  in  the  object  plane.  The  second  source  of  degradation 
is  an  additive  noise  represented  by  n(*  )  which  can  only  be  character¬ 
ized  in  statistical  terms.  Figure  2-1  represents  a  linear  imaging 
system  and  the  corresponding  block  diagram.  Generally,  the  response 
h(-)  in  the  image  space  varies  with  the  position  ( ^ ,  Tj )  of  the  input 
impulse  and  is  called  a  space-variant  point-spread  function  (SVPSF) 
in  an  optical  context.  If  h( * )  is  isoplanatic,  i.  e. ,  the  form  of  h(.  ) 
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remains  fixed  in  the  image  plane  for  all  (  £  ,  n  )  positions, 
then  the  system  is  said  to  be  spatially  invariant  and  h(*  )  is  called  a 
space-invariant  point-spread  function  (SlPSF).  In  this  case,  the 
impulse  response  is  a  function  of  two  variables  and  the  dimensionality 
of  the  system  reduces  considerably.  The  PSF  h(x,  y;~,T  )  can  then  be 
written  as  h(x-c,y-Tl)  and  the  superposition  integral  (2.1)  simplifies 
to  a  convolution. 

CD 

g(x,  y )  =  JJ  h(x-c  ,  y-D)f(?  ,  T'l)ds'din  +  n(x,  y).  (2.2) 

-00 


The  mathematical  representation  given  in  (2.2)  is  general  enough  to 
cover  many  situations  that  occur  in  coherent  and  incoherent  optical 
systems . 

Some  of  the  sources  of  degradation  include:  diffraction,  motion 
degradation;  defocusing  and  atmospheric  turbulence.  Diffraction 
in  an  optical  imaging  system  is  due  to  the  limited  aperture  size  and 
is  an  example  of  spatially  invariant  degradation.  The  blur  function 
for  a  system  with  a  circular  aperture  and  incoherent  illumination  is 
given  by  f2 3 3 


(2.  1) 


where  p  =  (x  +y  and  Jj  is  a  Bessel  function  of  the  first  kind, 
order  one.  Linear  uniform  motion  degradation,  defocusing  and 
atmospheric  turbulence  are  other  examples  of  space-invariant 
degradations  [4*]  ,  [23"!  .  In  some  cases,  such  as  defocusinn  and 


and  atmospheric  turbulence,  the  impulse  response  is  separable,  i.e., 
the  function  of  two  variables  can  be  written  as  a  product  of  two 
functions,  each  with  one  variable. 

The  assumption  of  space-invariance  is  not  valid  for  certain 
degradations.  Lens  aberrations  such  as  coma,  astigmatism, 
curvature-of-field,  motion  blur  where  objects  are  at  different  distances 
from  the  camera  and  image  plane  tilt  are  examples  of  space -variant 
systems  [26*]  .  By  an  appropriate  selection  of  coordinates,  some  of 
these  degradations  can  be  transformed  into  equivalent  space- 
invariant  [24l  ,  [25l  systems. 

The  assumption  of  additive  noise  is  broad  enough  to  encompass 
different  practical  situations.  Many  of  the  noise  sources  (e.g.  , 
stray  illumination,  circuit  noise,  roundoff  error)  may  be  individually 
modeled  as  additive  noise.  Nevertheless,  the  assumptions  of  linearity 
and  additive  noise  are  subject  to  criticism  because  they  are  valid  only 
over  a  certain  dynamic  range.  The  problem  is  that  g  is  not  directly 
available  for  processing.  Instead,  a  nonlinear  recording  of  g  on  a 
photographic  emulsion  is  usually  the  only  available  measurement.  It 
is  possible  to  measure  the  nonlinear  function  to  recover  g  over  a 
larger  dynamic  range,  but,  any  attempt  at  extending  this  range  must 
ultimately  be  frustrated  by  a  drastic  increase  in  the  noise  level  f4"l  . 
Also,  the  effect  of  film  grain  noise  is  far  from  being  additive. 

Huang  [27*1  has  shown  it  could  be  modeled  by  a  multiplicative  process, 

1  1 


and  more  general  signal-dependent  models  must  be  used  to  accurately 
describe  the  process  [28l  . 

After  specifying  assumptions  and  limitations,  the  next  step  is  to 
clarify  the  necessary  information.  The  model  assumes  that  a 
complete  knowledge  of  the  impulse  response  h  is  available.  This 
knowledge  can  be  obtained  analytically  [23"'  or  from  edges  or  points 
in  the  image  that  are  known  to  exist  in  the  object  [29l  .  As  far  as 
the  noise  is  concerned,  knowledge  of  the  second  order  statistical 
properties  is  required.  The  noise  is  not  necessarily  white,  but  this 
assumption  is  often  made. 

Each  restoration  scheme  given  in  the  succeeding  sections  and 
chapters  assumes  some  objective  intuitively  reasonable  criteria  of 
quality.  Inverse  filtering,  for  instance,  attempts  perfect  resolution 
without  regard  to  noise,  while  the  Wiener  filter  minimizes  the  mean 
square  error  without  regard  to  resolution.  Although  it  is  known  that 
the  human  observer  does  not  judge  images  according  to  mean  square 
error  [30],  it  has  been  found  that  reasonable  results  canbe  obtained  by 
its  use, especially  for  low  contrast  images.  Moreover,  mean  square 
error  leads  to  a  very  tractable  mathematical  structure,  the  regres¬ 
sion  model,  which  has  been  considerably  explored  in  mathematical 
statistics. 

2.  2  Inverse  Filtering 

The  idea  of  inverse  filtering  is  very  simple.  Taking  a  Fourier 
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transform  of  both  sides  of  the  convolution  expression  (2.2)  gives 

G(u,  v)  =  H(u,v)F(u,v)  +  N(u,  v)  (2.4) 

where  u  and  v  represent  the  x  and  y  spatial  frequencies  and  the 
upper-case  letters  represent  Fourier  transforms  of  the  functions 
denoted  by  the  corresponding  lower-case  letters.  If  the  transfer 
function  H(u,  v)  does  not  vanish  at  any  point  (u,  v),  the  inverse  filter 
R(u,v)  is  defined  as 

R(u,  v)  =  1/H(u,v)  .  (2.5) 

The  restored  image  in  the  Fourier  domain  is  obtained  by  multiplying 
both  sides  of  (2.4)  by  R  and  taking  an  inverse  Fourier  transform. 
Thus 


A 

F(u,  v) 

=  F(u,  v)  +  N(u,  v)/H(u,  v) 

(2.6) 

A 

f(x,  y) 

=  f(x,  y)  +  jVn/H) 

(2.7) 

where  represents  the  inverse  Fourier  transform.  Marechal 
etal.  [3],  Tsujiuchi  [5*1  ,  Harris  [6*1,  McGlamery  f7]  ,  and  Mueller 
and  Reynolds  [8*1  h  -e  applied  this  method  with  minor  modifications. 
These  modifications  aim  at  two  short-comings  of  the  inverse  filter. 
First,  since  H  in  most  practical  cases  decreases  rapidly  for  large 
values  of  u  and  v,  and  the  noise  has  a  fairly  flat  spectral  density 
function,  this  filter  amplifies  the  high  frequency  noise.  Second,  in 


many  cases  H  has  zeros  at  spatial  frequencies  within  the  range  of 
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interest,  and  since  division  by  zero  is  meaningless,  the  method 


breaks  down.  Therefore,  a  perfect  restoration  is  impossible  even  in 
the  noiseless  case.  All  the  modifications  suggested  to  overcome 
these  short-comings  are  intuitive  and  ad  hoc.  For  example,  Harris 
[6"1  suggests  that  the  right  hand  side  of  (2.5)  should  be  multiplied  by 
a  function  H^lu,  v)  before  taking  the  inverse  Fourier  transform.  11^ 
is  chosen  to  be  zero  wherever  the  power  spectral  density  of  the  signal 
is  dominated  by  noise  and  II ^ /H  is  chosen  to  be  finite  at  the  zeros  of 
H.  With  this  method,  the  problem  of  infinities  is  removed,  but  the 

A 

restored  image  is  a  blurred  version  of  f  derived  in  (2.  6)  through  the 
filter  Hq.  The  requirements  imposed  on  are  intuitively  reason¬ 
able.  However,  there  are  an  infinite  number  of  functions  satisfying 
these  requirements  and  one  of  them  must  be  chosen  arbitrarily. 

It  has  been  shown  that  with  the  modifications  given  in  mrsi-tsi. 
and  in  spite  of  the  limitations  and  assumptions,  inverse  filtering  can 
yield  reasonably  good  restoration  in  many  cases  where  noise  is  not 
dominant. 

2.  3  Wiener  Filtering 

The  basic  idea  of  Wiener  filtering  is  to  minimize  the  mean 
square  error  (MSE)  between  the  original  and  restored  images.  MSE 
is  an  objective  criteria  for  which  the  optimum  restoration  can  be 
rigorously  computed.  The  following  is  a  '  **ief  outline  of  the  mathe¬ 
matical  derivation  of  the  Wiener  filter. 
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Assuming  the  object  function  f  and  the  noise  n  be  sample  func- 
tions  of  different  random  processes,  the  minimum  MSE  estimate  of 
f  is  f  such  that  for  each  (x,y)  of  the  object  plane  the  error  e  given  by 

e  =  E[(f-f)21  (2.8) 

is  minimized,  where  E  denotes  expectation  over  the  processes  f  and 
n.  It  can  be  proven  that 

f  =  E(f|g)  (2.9) 

in  general. 

Although  (2.9)  appears  very  simple,  the  computation  of  the 
conditional  expectation  is  very  difficult  for  arbitrary  processes  f 
and  n.  The  problem  becomes  much  simpler  if  we  assume  these 
processes  to  be  stationary,  independent,  zero-mean  Gaussian 
processes.  In  this  case,  it  can  be  shown  that  the  optimum  nonlinear 
or  linear  estimate  of  (2.9)  can  be  obtained  by  linear  spatially  in¬ 
variant  filtering  of  g.  This  minimization  filter  is  easier  to  des¬ 
cribe  in  the  spatial  frequency  domain,  and  can  be  shown  to  be  [4*1 

R  =  H*/(HH*  +  cp  /cp,)  (2.10) 

n  i 

where  *  denotes  complex  conjugate  and  Cp^  and  are  power  spectral 

densities  of  the  processes  f  and  n.  With  cp^  =  0,  the  Wiener  filter 

becomes  the  inverse  filter  and  when  H  =  0,  R  =  0  (unlike  R  =  •  with 

the  inverse  filter).  In  order  to  determine  R,  both  H  and  the  noise-to- 
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signal  ratio  cp  /tp,  must  be  known.  In  most  cases  33  /cp,  is  assumed 
^  n  f  n  1 


to  be  constant  over  the  frequency  range  of  interest. 

For  a  better  fidelity  criterion  than  MSE,  various  modifications  to 
the  Wiener  filter  have  been  proposed.  Hunt  [14*1  has  proposed  a 
constrained  least-square  filter  and  Stockham  and  Cole  4  5  '  have 
suggested  a  filter  which  is  the  geometrical  mean  between  the  inverse 
filter  and  the  Wiener  filter.  These  filters  have  been  implemented 
with  very  successful  results. 
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Chapter  3 


DISCRE  TIZATION  OF  THE  CONTINUOUS  MODEL 

In  the  processing  of  images  by  digital  computer,  the  continuous 
model  of  Eq.  (2.  1  )  must  be  discretized.  In  digital  image  processing, 
the  information  is  necessarily  finite  and  discrete  in  both  amplitude 
and  spatial  position.  Therefore,  the  continuous  image  field,  and  in 
most  cases  the  impulse  response,  must  be  transformed  into  arrays 
of  numbers.  Generally,  this  transformation  produces  some  error, 
i.e.,  the  inverse  transform  of  these  arrays  of  numbers  is  not 
exactly  the  original  image  field. 

Representing  the  continuous  function  by  an  array  of  samples, 
known  as  the  pulse  approximation  method,iS  the  simplest  and  the  most 
common  technique  in  image  discretization.  However,  the  accuracy 
of  this  technique  is  suspect  in  digital  image  modeling.  Using 
numerical  analysis  methods,  such  as  quadrature  formulae,  leads 
to  a  more  accurate  model.  Spline  functions,  because  of  their  highly 
desirable  interpolating  and  approximating  characteristics,  are 
suggested  as  a  potential  alternative  to  the  conventional  pulse  approxi¬ 
mation  method. 

In  this  chapter,  the  problems  of  image  sampling  and  quadrature 
formulae  are  analyzed.  It  is  shown  that  spline  functions  are  superior 
to  both  pulse  approximation  techniques  and  polynomials  in  discrete 
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representation  of  a  continuous  function  and  numerical  solution  of 
integral  equations.  Some  experimental  results  are  given  in  the  last 
section  of  this  chapter. 

3.  1  Pulse  Approximation  Method 

The  idea  of  the  pulse  approximation  method  is  to  represent  a 
function  f(x)  by  an  array  of  its  sampled  values  taken  on  a  countable  set 
of  points  on  the  x  axis.  Clearly,  if  the  sample  points  are  close 
enough,  the  sampled  data  are  an  accurate  representation  of  the 
original  picture.  Thus  the  function  f  can  be  reconstructed  with  suf¬ 
ficient  accuracy  by  simple  interpolation.  Assuming  Ax  to  be  the 
distance  between  two  subsequent  points  in  a  uniform  sampling,  the 
sampled  function  fg(x)  is  obtained  by  multiplying  the  original  function 
by  summation  of  6  functions  as  expressed  by 

OP 

f  (x)  =  E  f(iAx)6(x-i Ax)  .  (3.1) 

i=-® 

Taking  a  Fourier  transform  of  (3.  1),  the  spectrum  of  the  sampled 
function  is  given  by  [23*] 


F  (u)  =  E  F(u  -  —  ) 
s  .* — *  Ax 

i  =  -» 


(3.2] 


where  F  is  the  Fourier  transform  of  f  and  u  is  spatial  frequency. 


Equation  (3.2)  indicates  that  the  spectrum  of  the  original  function  is 
infinitely  repeated  with  a  distance  of  .  Assuming  the  function  f 


to  be  bandlimited,  its  spectrum  F  is  non-zero  over  only  a  finite 
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interval  R  in  the  frequency  domain.  If  Ax  is  sufficiently  small,  then 


the  separation  — —  is  large  enough  to  assure  that  adjacent  spectra  do 
not  overlap.  If  2B^  represents  the  width  of  the  rectangle  that 
completely  encloses  the  interval  R,  then  non-overlapping  is  assured 
if 


Ax 


(3.  3) 


x 

Physically,  this  means  that  the  function  f  must  be  sampled  at  a  rate 
at  least  twice  its  highest  frequency  component  or  one-half  the  period 
of  the  finest  detail  within  the  function.  If  equality  holds  in  l'q.  (3.  31, 
the  function  is  said  to  be  sampled  at  its  Nyquist  rate.  If  Ax  is 
smaller  or  larger  than  this  threshold,  the  function  is  oversampled  or 
under  sampled.  With  condition  (3.3)  satisfied,  the  exact  reconstruc¬ 
tion  of  the  original  function  can  be  achieved  by  filtering  the  sampled 
data  with  an  appropriate  filter,  for  example  a  filter  with  a  rectangular 
transfer  function  of  width  2B^.  In  the  spatial  domain,  the  reconstruc¬ 
tion  operation  in  the  spatial  domain  is 


f(x)  =  f(2B~  )sinC^2BxfX_2:B_  ^  (3-4' 

i =-*  x  x 

for  a  rectangular  filter.  Equation  (3,4),  known  as  the  Whittaker- 
Shannon  sampling  theorem,  indicates  that  the  function  is  reconstructed 
exactly  by  an  infinite  sum  of  weighted  sine  functions  injected  at  each 
sample  point. 
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3.  2  Quadrature  Formulae 


The  problem  of  image  degradation,  as  stated  in  Eq.  (2,  1),  is 

represented  by  an  integral  equation.  In  practical  situations,  the 

limits  on  this  definite  ntegral  equation  are  not  infinite.  First,  the 

degradation  function  h  usually  vanishes  (or  almost  vanishes)  beyond 

some  point,  and  consequently,  h  is  non-zero  over  a  finite  interval. 

Second,  only  a  finite  size  of  the  object  is  of  particular  interest  for 

restoration.  With  these  considerations,  the  one -dimensional 

version  of  Eq.  (2.  1)  is  a  definite  integral 

b 

g(x)  =  J  h(x,  r)f(“  )dr  (3.5) 

a 

over  a  finite  interval  [a,  bl  . 

To  implement  this  continuous  integral  by  a  digital  computer,  a 
numerical  technique,  called  a  quadrature  formulae  (q.  f.  )  must  be 
employed.  A  q.  f.  is  an  approximation  to  a  definite  integral  by  a 
linear  combination  of  values  of  the  integrand,  and  perhaps  also  of 
some  of  its  derivatives,  at  certain  points  of  the  interval  of  integra¬ 
tion  called  the  nodes  of  the  q.  f.  [31 1.  A  discrete  version 


g(x. ) 


c  h(x., 

i.)  1 


r.)f(c.) 
3  3 


(3.  6) 


of  Eq.  (3.  5)  can  be  obtained  by  applying  a  q.  f.  Using  vector  space 
notation,  the  above  equation  simplifies  to 


£  =  Hf 


(3.7) 
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where  II  is  an  M  x  N  matrix  with  elements  h  -  e  h(x  ,  '  I. 

—  u  1.1  1  1 

Assuming  the  coefficients  e  =  1  is  equivalent  to  the  pvilse  approxima¬ 
te 

tion  method.  For  a  given  number  of  samples,  a  good  choice  of  q.  f. 
can  result  in  an  accurate  vector  space  model.  Moreover,  the  quadra¬ 
ture  coefficients  can  affect  the  stability  of  the  model  and  decrease 
the  condition  number. 

The  general  form  of  a  q.  f.  ,  when  the  derivatives  are  not  avail¬ 
able,  is  given  by 

,b  XT'* 

j  f(x)dx  =  c.f(x.)+  Rf  (3.  hi 

a  i=l  ’  ’ 

where  c.  and  x.  are  coefficients  and  nodes  of  the  q.  f.  ,  respe  ctively, 
r he  term  Rf  is  a  functional  which  for  any  given  function  ft*  )  equals 
the  difference  between  the  exact  value  of  the  integral  and  its  approxi¬ 
mation.  For  a  given  q.  f.  ,  Rf  depends  on  the  integrand  and  may  vanish 
for  some  specific  class  of  functions.  Therefore,  the  objective  is  to 
minimize  the  upper  bound  of  Rf  as  well  as  to  enlarge  the  class  of 
functions  which  result  in  zero  value  for  Rf.  If  the  nodes  of  the  q.f. 
are  pro -assigned,  the  only  available  parameters  to  be  treated  are 
the  coefficients.  Examples  of  this  type  are  Newton-Cotes  and  best 
q.f.  in  the  sense  of  Sard  [31 1  .  If  the  nodes  are  free,  the  best 
location  of  the  nodes,  in  a  certain  sense,  can  be  determined,  and  the 
q.f.  is  called  optimal.  Examples  of  the  optimal  type  are  Gauss- 

Legendre  and  optimal  q.f.  in  the  sense  of  Sard.  Since  in  most  cases, 
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particularly  in  image  processing,  the  location  of  the  nodes  are  pre¬ 
assigned,  only  fixed  node  q.f.  are  considered  here.  Newton-Cotes 
q.  f.  is  briefly  studied  in  this  section,  best  q.f.  in  the  sense  of  Sard 
in  section  3.5  and  the  experimental  results  are  compared  in  section 

3.6. 

The  basic  idea  of  Newton-Cotes  q.f.  is  to  interpolate  the 
sampled  data  by  Lagrange  method  and  then  integrate  it  [32"1.  Clearly 
the  remainder  of  the  integral  has  the  property  that  RF  =  0  if  fetr 

n-  i 

where  n  is  the  entire  class  of  polynomials  of  degree  less  Lhan  or 
equal  to  n-1.  This  property  may  be  used  to  determine  the  coeffi¬ 
cients  c,  .....  c  .  A  linear  system  of  equations  can  be  obtained  by 
In 

n-1 

assuming  Rf  -  0  when  f(x)  =  1 ,  x,  .  .  . ,  x  in  equation  (3.  8).  The 
coefficients  are  the  solution  of  this  linear  system  of  equations. 

3.  3  Spline  Functions 

Spline  functions  are  a  class  of  piecewise  polynomial  functions 

satisfying  continuity  properties  only  slightly  less  strigent  than  those 

of  polynomials,  and  thus  they  are  a  natural  generalization  of 

polynomials  [3 3 "1  .  Given  a  strictly  increasing  sequence  of  real 

numbers  x  ,  x  ,  . . .  ,  x  ,  a  spline  function  S(x)  of  degree  m  with  the 
1  2  n 

knots  x, ,  x  ,  . . . ,  x  is  a  function  having  the  following  two  properties; 

1  2  n 

1)  In  each  interval  (x^,x. S(x)  is  given  by  some  polynomial 
of  degree  m  or  less. 

2)  S(x)  and  its  derivatives  of  order  1 , 2,  . . . ,  m-1  are 
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continuous  everywhere. 


When  m  =  0,  Condition  2  is  not  operative,  and  a  spline  function  of 
degree  zero  is  a  step  function.  A  spline  function  of  degree  one  is  a 
polygon. 

In  general,  the  polynomials  representing  S(x)  in  adjacent  inter¬ 
vals  (x.,  x.  ,  )  and  (x  ,  ,  x  _)  are  different,  although  this  is  not  a 

i  t+1  i+l  i+2 

requirement.  S/x)  might  be  represented  by  a  single  polynomial  on 
the  entire  real  line.  In  other  words,  all  the  polynomials  of  degree 
m  or  less  are  included  in  the  class  of  spline  functions  satisfying  the 
above  properties.  Spline  functions  can  equally  be  defined  as  the 
following: 

1 )  For  m  >  0,  a  spline  function  of  degree  m  is  a  function  in  the 

class  of  m-1  times  differentiable  functions  (C™  )  whose 

derivative  is  a  step  function. 

th 

2)  A  spline  function  of  degree  m  is  any  m  order  indefinite 
integral  of  a  step  function. 

Polynomials, because  of  their  simple  mathematical  properties, 
have  been  widely  used  for  interpolation  and  approximation.  However, 
a  polynomial  fitted  to  a  fairly  large  number  of  data  points  has  numer¬ 
ous  and  severe  undulations.  There  is  now  considerable  evidence  that 
spline  functions  in  many  situations  are  more  adaptable  approximating 
functions  than  polynomials  with  a  comparable  number  of  parameters. 


Moreover,  they  havebeen  showntobe  the  solution  of  some  optimization 
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problems  f34"!  ,  [35l  ,  [3(>1  . 

I'he  basis  for  tho  class  of  spline  functions  of  degree  m  having 

knots  . . x  is  given  by 

14m 


fl, 


m  m  ,m 

X  ,  (X-Xj  1  ,  (x-x^  , 


.m-, 

(x-x  )  ^ 

n  + 


(3.9) 


where 


m 


(x-x.) 

i  + 


0 


if  x  £  x. 

i 


(3.  10) 


(m 

(x-x. )  if  x  >  x. 
i  i 

Using  this  basis  for  interpolation  and  approximation  turns  out  to  be 

unstable  in  practice,  since  the  matrix  of  the  system  is  very  badly 

conditioned  unless  m  and  n  are  both  small  [37l  .  The  numerical 

instability  is  related  to  the  mathematical  properties  of  the  truncated 

power  functions.  This  difficulty  can  be  overcome  by  adopting  another 

basis  for  the  class  of  spline  functions.  The  most  desirable  basis 

consists  of  splines  with  finite  support  containing  a  minimum  number 

of  knots.  This  basis,  called  B-splines,  has  minimal  support  for  a 

given  degree  and  has  been  studied  by  Curry  and  Schoenberg  [38*1  . 

A  B  -spline  M.  (x)  of  degree  m  with  knots  x  ,  x.  x.  ,  is 

i  i  i  +  l  i+m  +  1 

given  by 


i+m+1  (x-x.) 

M.W  =  (m+1)  £  -JjJJ 


(3.  11) 
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where 


i+rn+1 

u>(x^  =  fl)rn+  IT  (x4-*j)-  (3.12) 

j=i 

If  the  knots  are  uniformly  spaced,  B-splines  have  the  following 
properties : 

1)  The  functions  are  shift-invariant,  i.e. 

M.'x)  =  M.  ,  (x-x  ) 

i  i  -k  k 

2)  The  functions  are  strictly  positive. 

3)  The  convolution  of  two  B-splines  of  degrees  m  and  n 
yields  another  B-spline  of  degree  mtn+1. 

4)  The  functions  have  limited  support  and  thus  are  a  local 
basis . 

These  interesting  properties  can  be  fully  exploited  in  image 
processing.  Property  1  results  in  circulant  matrices  which  can  be 
inverted  by  a  Fourier  transform.  The  second  property  is  useful 
because  image  intensity  is  always  non-negative.  Property  3  is 
analogous  to  convolution  in  space -invariant  degradation  which  will  be 
the  subject  of  section  4.  1.  Finally,  the  fourth  property  produces 
banded  matrices  which  are  very  efficient  in  computation. 

Figures  3-1  and  3-2  are  plots  of  B  -  splines  of  degrees  1  ,  2  and  3 

centered  at  the  origin  of  the  coordinate  system.  A  B-spline  of  degree 

zero  is  a  rectangle  function  and  a  B-spline  of  degree  one  is  a  riangle 
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function.  Using  Eqs.  (3.  11)  and  '3.  12),  B  -splines  of  degrees  2  and  3 


can  be  expressed  as 


B2(x)  =  3 


(x-f  1 .  5)*^  (x  +  .  5)  +  (x-  .  5 )*f 


B  3  (x)  =  4 


[" (x+2  ) ^  _  (x+1  )  +  ^  (x)  +  fx-  1  )1  ^  fx-  2  )1 

~~ ZT~  6  4  '  6  24 


- 

] 


13) 


(3. 14) 


Spline  functions  of  degree  three,  called  cubic  splines,  in  many 
situations  have  more  desirable  properties  than  other  splines,  and 
therefore  they  are  widely  used  for  approximation  and  interpolation. 

3 . 4  Error  Analysis 

Interpolation  and  approximation  by  spline  functions  is  generally 
not  without  error.  For  a  given  function,  the  error  depends  on  the 
function,  the  degree  of  spline,  and  the  number  of  placement  of  the 
knots.  An  analysis  of  the  error  is  helpful  in  choosing  the  proper 
spline  function.  The  following  is  a  brief  error  analysis. 

Let  PC™’  (a,b)  be  the  set  of  all  real- valued  functions  fix)  such 

that : 


1) 


2) 


3) 


f(x)  is  m-1  times  continuously  differentiable  on  the  open 
interval  (a,  b). 

There  exist  a  sequence  of  knot  s  a  =  x„  <  x  <  x  <  .  .  .  <  x  < 

0  12  n 

x  ,  =  b  such  t  hat  on  each  open  interval  fx. ,  x.  ,),  0  i  i  4n, 
n  f 1  ii  +  l 

f  is  m  times  continuously  differentiable. 

The  Lp-norm  of  derivative  is  finite,  i .  e.  , 


n  i  4 1  n  1  /P 

Dmf  p  =  j  |Dmf(x)fdxj 


<  »  .  (3.15) 


i  =  0  x. 

When  p  =  «,  the  third  requirement  becomes 


I!  DmfH  =  max  sup  |nmf(x)l<«.  (3.16) 

oo 

O^i^n  xefx. , x  .  ) 
i  i+l 

With  the  above  definitions,  Shultz  [39*1  has  derived  error  bounds 
for  different  functions  interpolated  by  cubic  splines.  Assuming  f(x) 
to  be  the  original  function  and  S(x)  to  corresponding  cubic  soline 
function,  the  error  bound  is  given  as  follows. 

2  2 

If  fcPC  ’  (a.b),  then 

I'f-Sl'  *  2TT'2h2|lD2f'2  .  (3.17a) 

2  0D 

If  fePC  ’  (a.b)  then 

Ilf-S'l  s|h2||D2f|!  .  (3.17b) 

00  3  00 

4  2 

If  fcPC  ’  (a.b),  then 

l|f-Sll  s4n‘4h4llD4f|L  .  (3.17c) 

If  fePC4,"(a,  b),  then 

in-sit  *  3§J  I.4!! o4f"..  n.lTd) 

Inequalities  (3.  I7a-d)  indicate  that  the  error  bound  is  a  mono¬ 
tonic  increasing  function  of  sampling  interval  h  and  the  norm  of  the 
m-th  derivative  of  the  original  function.  Therefore,  one  way  to 
reduce  the  error  is  to  sample  the  function  at  a  higher  rate.  When 
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the  function  f  belongs  to  more  than  one  of  the  above  categories,  the 
error  bound  is  the  smallest  one.  If  f  is  a  polynomial  of  degree  less 
than  or  equal  to  3,  then 

"  D4f  =  ilD4f"  =  0  , 

2  00 

therefore* 

i'f-si'  =  I'f-s!1  =  o 

2  » 

and 


f  (x )  =  S(x) 


which  means  cubic  splines  exactly  interpolate  the  polynomials  of 
degree  3.  As  another  example,  let  f(x)  be  a  sine  function  with 
frequency  u,  then 

2  2 
D  f  =  -(2ttu)  sin  2ttux 


and 


Therefore 


and 


D4f 


(2hu)  sin  2itux  . 


!'D2f!'  =  (2nu)2 

’  OB 

"D4fl!  =  (2ttu)4  . 

os 


Substituting  the  above  norms  in  <3.  I7b,d) 
||f-S|l  s  (2nhu)2 

and 
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Hu-  error  bound  is  the  minimum  of  the  above  t\v<.  limits.  For  a 


given  error  bound,  h  is  proportional  to  the-  inverse  of  u  which  is 

similar  to  the  sampling  theory  studied  in  section  3 .  1. 

3.  5  Monospline  and  Best  Quadrature  Formulae  in  i  he  Sense  of  Sard 

Monosplines  are  a  class  of  functions  defined  as  [3l] 

m 


K(x)  = 


-  S 


(x) 


where  S 


m  '  m  -  1  ,  n 
(x)  is  a  spline  function  of  degree  m-1  with  n  pre- 


(3.  181 


m  - 1 ,  n 

assigned  knots  a  <  x,  <  x_<.  .  .  <  x  <  b  and  n  >  m  given  by 

1  2  n 

m-1 


S  .  fx) 
m- 1 ,  n 


m-1 

=  £ 

j=0  3  i=l 


n  c.(x-x. 


i  4 


(m-  I  1 ! 


(3.19) 


K(x)  which  consists  of  a  polynomial  of  degree  m  and  a  spline  of 
degree  m-lis  called  a  monospline  of  degree  m  with  n  nodes.  Using 
K  (x)  as  a  kernel 


)dx 


J  f(x)K^m  (x)dx  =  J  f(x)dx  -  j*  c.f(x)6(x-x. 
a  a  i  =  l  'a 

b  n 

=  j  f(x)dx  -  c.f(x.)  . 

a  i  =  1 

Integrating  the  left  hand  side  of  (3.20)  m  times  (by  parts)  gives 

J  f(x)K(m)(x)dx  =  £  (-  1  )3f^3 )(x)K^m" 1  ~3\x) 
a  j  =  0 


(3.20) 


f,m,(x)K(xldx  .  (3.21^ 


L. 


Assuming 


Kfj)(a)  = 


K(^(b)  = 


0  for  j  =  0,1, 


m  - 1 


(3.22) 


and  substituting  (3.21)  in  (3.20)  gives 

b  n  .b  . 

J  f(x)dx  =  ^  c.f(x.)  +  1  f  m  (x)K(x)dx.  (3.23) 

a  i=l  1  a 

Therefore  the  remainder  of  the  integral,  or  the  error  of  q.f.  ,  is 
given  by 


Rf 


f^m^(x)K(x)dx 


a 

The  upper  bound  of  Rf  can  be  expressed  as 


(3.24) 


|  Rf  |  =  |  ^  f(m)(x)K(x)  |  s  !iK!!2!lf(m)l'2  (3.25) 

a 

where  ||*H_  denotes  L  -norm  of  the  function.  If  Rf  =  0  when  f  is  a 
£•  £* 

polynomial  of  degree  less  than  or  equal  to  m  1  and  K(x)  has  the  least 

square  deviation  (minimum  norm)  among  all  kernels  of  the  form 

(3.  18),  then  the  q.f.  is  called  the  best  in  the  sense  of  Sard.  For  a 

given  function  f,  the  minimum  norm  of  K  generates  the  minimum 

upper  bound  of  Rf.  Assumption  (3. 22)  which  leads  to  Eq.  (3,  23) 

satisfies  the  first  requirement.  Schoenberg  [40"]  ,  [41 1  has  shown 

that  there  exists  a  unique  monospline 

2m 


H(x) 


x _ 

(2m)! 


-  S  (x) 

2m- 1 ,  n 


(3.26) 


of  degree  2m  in  which  the  kernel  K(x)  of  Sard's  best  q.f.  ,  in  terms  of 
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H(x),  is  given  by 


K(x)  =  H(m)(x)  .  (3.27) 

H(x)  must  satisfy  the  following  conditions 


H(x.)  =  0  i  =  1,2,,.. ,n  (3.28a) 

i 

H^rn+j^(a)  =  0  j  =  0 .  1 . m-1  (3. 2  8b) 

H(m+j\b)  =  0  j  =  0,  1 . m-1  (3.28c) 


! ,  C, 

The  value  J  of  minimum  derivation  ,'K  'n  terms  of  H(x),  is 
determined  by  the  relation 


Assuming  a 
the  normalizing 
be  written  as 


b 

m  (* 

( - 1 )  J  H(x)dx  . 


(3.29) 


a  a 

=  -I  ,  b  =  1  (this  assumption  can  always  be  made  by 
S  =  — —  ),  and  applying  condition  (3.28b),  H(x)  can 


I!(x) 


(x+1) 

(2m) 


2m  m-1  _n  c.(x-x 

-  Sv'  l  ■** 


2m  -  1 

) 

4 

-1  )' 


(3.  30) 


Conditions  (3,28a)  and  (3.28c)  generate  a  system  of  min  linear 


equations  with  m+n  unknowns  of  the  form 


m-1  . 

(2m  II  '  S  aiX'i 

J  =  0  • 


c.(x.  _x. ' 

.  2—1 _ J_ 


,2m  -  1 


(2m-  1 ) 1 


=  0 


i  =  1  ,  .  .  .  ,  m 


(3.  31a) 
(3.  ’.lb) 
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The  coefficients  c.  and  a  and  the  L  -norm  of  K  can  be  obtained 

.1  j  2 

by  solving  the  above  linear  system  of  equations.  If  m  =  1,  which 
corresponds  to  the  pulse  approximation  method,  the  system  can 
easily  be  solved.  In  this  case  the  equations  are  given  by 


(x.  t  l  ) 

i 


2  '  a0 


n 

-  c.(x  -x  )  =0 

U  -1  1 


1,2 . n 


II 

S'.  ■ 


2  . 


The  solution  to  the  above  equation  is  obtained  as 

2 


0 


1 1  t  x  j  ) 


2 +xi ,x2 


x.  .  -x. 

J  + 1  J  !  1 


i  -  2 ,  1 . n  -  1 


2  -x  . - x 
n  -  1  n 


'  3.  32a) 


M.  32b) 


i 3 .  33a  ) 

i  3.  3  3b) 

( 3.  33c) 

l 3.  33d) 


When  the  sample  points  arc  equidistant,  the  location  of  nodes  and 
coefficients  are  derived  as  follows 


x. 

.1 


-1  + 


2i-l 


i  =  1,2 . n  (3.  34a  ) 


c . 
.1 


2 

n 


1 


i  =  1 , 2,  .... n  i  3 .  3 4 b ) 


f  3.  34c  ) 


Substituting  13.  34a -c)  in  (  3.  30)  produces 


H(x)  = 


(x+1)  1 


7  -  -  (x+x.) 

2„2  "j  =  l  J  + 


(3. 35) 


The  value  J  of  mini  mum  deviation  can  be  obtained  as 

1  1 
J  =  J  (K(x))  dx  =  (l-)mJ  H|x)dx 


■  1 


1 


=  -f‘  d*  +  1 

-1  L 


\  *  ->  n  1 

- j  f  dx  +  —  f  (x-x. 

2n2  > 


-x.)  dx 

+ 


8  1  1  *7*  1  2j  2 

—  f  —  +  —  y.  (2  +  —  ) 

b  £■  n  n  n 

n 


3n 


2  * 


Substituting  the  norm  of  K(x)  in  (3.25),  the  upper  bound  of  the  error 


is 


R(s  i  /fllf'l 


(3.  3 b) 


Therefore,  as  was  expected,  the  upper  bound  of  error  is  inversely 
proportional  to  the  number  of  sample  points  and  approaches  zero  as 
n  increases.  Moreover,  Rf  =  0  for  constant  functions  f(x)  =  c 
regardless  of  the  number  of  sample  points. 

3.  6  Experimental  Results 

To  show  the  improvements  that  can  be  made  by  using  mono¬ 
splines,  this  section  is  devoted  to  applying  Sard's  best  q.f.  to  a 
variety  of  functions  and  comparing  the  results  with  pulse  approxima¬ 
tion  method  and  Newton-Cotes  q.f. 


3  5 


In  applying  Sard's  best  q.  f.  ,  one  is  faced  with  the  task  of 


selecting  the  parameter  m.  For  a  given  number  of  sample  points, 

the  deviation  of  K(x)  decreases  as  m  increases.  Figure  3-3  is  an 

illustration  of  this  property  for  8  uniformly  spaced  nodes.  There- 

f  m ) 

fore,  one  may  assume  m  to  be  the  highest  order  where  f  is 
continuous.  On  the  other  hand,  as  m  increases,  different  problems 
will  arise.  First,  the  system  of  equations  (3.  3  1  a,  b)  tends  to  become 
unstable  for  large  m.  Second,  in  some  cases,  the  norm  of  the  rr.-th 
derivative  of  the  integrand  increases  rapidly  as  m  grows,  and  a 
smaller  choice  of  m  would  result  in  a  smaller  error.  To  study  the 
effects  of  different  values  of  m  on  the  error  and  also  to  compare 
Sard's  method  to  the  Newton-Cotes  and  pulse  approximation  methods, 
several  experiments  have  been  performed.  Figure  3-4a  demonstrates 
the  error  as  a  function  of  frequency  for  a  sine  function.  Case  m  =1 
coincides  with  the  pulse  approximation  method,  and  m  =  8  is  equi¬ 
valent  to  the  Newton-Cotes  q.  f.  Figure  3-4b  is  a  plot  of  the  theoreti¬ 
cal  error  bound  for  a  sine  function.  In  Figure  3-5a  the  integrand  is 
a  polynomial  of  degree  8.  Variable  j  is  a  measure  of  how  fast  the 
polynomial  oscillates  in  the  interval  of  integration;  j  is  almost  equal 
to  the  number  of  roots  in  the  interval  (-1,  1)  minus  one.  In  other 
words,  the  larger  j  becomes,  the  harder  it  is  to  approximate  the 
function  since  it  is  subject  to  more  fluctuation.  This  roughly 


corresponds  to  the  frequency  in  Fig.  3-4.  Figure  i-Sh  shows  the 


Figure  3-3 


Figure  3-5.  Quadrature  error  for  polynomials. 


I’rror  for  polynomials  of  different  degrees.  In  this  plot  all  the  roots 
arc  between  -1  and  1. 

Both  theory  and  experience  indicate  that  the  choice  of  m  greatly 
depends  on  the  frequency  content  of  the  integrand  f.  For  the  class  of 
rapidly  varying  functions,  a  smaller  m  is  advised,  but  for  the  class 
of  slowly  varying  functions,  large  values  of  m  give  better  results. 
Since  it  is  assumed  that  f  is  sampled  faster  or  equal  to  the  Nyquist 
rate,  the  curves  of  Fig.  3-4  are  not  studied  for  frequencies  above 
two  cycles.  Figures  3-4  and  3-5  show  that  the  curves  cross  each 
other  and  that  a  tradeoff  exists  between  the  frequency  content  of  the 
integrand  and  degree  of  monospline.  Considering  this  fact  and 
taking  into  account  the  set  of  examples,  the  cubic  monospline  produces 
less  error  overall  and  thus  the  optimal  value  for  m  is  three.  Of 
course,  for  other  cases  where  the  function  f  is  highly  oversampled, 
large  values  of  m  may  be  recommended,  while  on  the  other  hand, 
when  the  function  is  sampled  far  below  the  Nyouist  rate,  the  pulse 
approximation  is  preferred  to  the  other  techniques. 

In  section  4.  4,  Sard's  best  q.  f.  has  been  used  in  the  simulation 
of  images  degraded  by  astigmatism  and  curvature  of  the  field.  This 
q.f.  results  in  a  more  accurate  model  with  the  reduction  of  simula¬ 
tion  artifacts.  Moreover,  this  technique  has  decreased  the  condition 
number  of  the  blur  matrix  and  consequently  produced  a  more  stable 
model  [24 1  ,  [42]  . 
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Chapter  4 


RESTORATION  OF  NOISELESS  IMAGES 


The  restoration  of  noise-free  images  is  presented  in  this 
chapter.  The  convolutional  property  of  B -splines  is  used  for  the 
restoration  of  space -invariant  degradations.  It  is  shown  that  repre¬ 
senting  the  object  and  point-spread  function  by  B-splines  leads  to  a 
more  accurate  reconstruction  of  the  original  object  than  the  conven¬ 
tional  method. 


The  singularity  of  most  imaging  systems  due  to  the  irreversible 
loss  of  original  object  information  is  a  major  problem  in  image 
restoration.  A  minimum  norm  principle  leading  to  pseudo-inversion 
is  used  to  overcome  this  difficulty.  This  technique  is  applicable  to 
space -variant  degradations,  underdetermined  models  and  over¬ 
determined  models.  The  space -variant  point-spread  functions  that 
describe  imaging  in  the  presence  of  astigmatism  and  curvature  of 
field  are  derived  and  coordinate  transformations  are  applied  to  reduc 
the  dimensionality.  The  singular -value -decomposition  is  used  for 
solution  of  the  simplified  equations. 


plication  of  B-splines  to  Fpace-Invariant  Degradations 


As  discussed  in  Chapter  2,  the  deterministic  part  of  a  degraded 


image  in  a  space-invariant  imaging  system  is  described  by  a 


convolution  integral.  Using  B-splines  as  a  basis  in  uniform 


4 


sampling,  the  object  f(xl  and  point-spread  function  h(x)  can  be 


represented  in  the  forms 


f  (x )  =  Y'  f.B  (x-x.) 

*—i  i  m  1 


'4.  1  ) 


h(x)  =  5  1  h.B  (x-x.) 

J  n  J 


(4.2) 


where  B  (x)  and  B  (x)  are  B-splines  of  degrees  m  and  n  centered  at 
m  n 

the  origin,  and  f.  and  h.  are  interpolation  coefficients.  Substituting 
(4.  1,2)  in  the  convolution  integral,  the  image  is 

00 

g(x)  =  j  h(x-?,)f(?)d" 


=  Vrf.h.B  (x-x.  )*B  (x-x.) 

i  j  m  i  n  i 


(4.  3) 


i  .1 


Exploiting  the  convolutional  property  of  B-splines 


B  (x-x.b;:B  (x-x.)  =  B  ,(x-x  -x.)  (4.4) 

m  i  n  j  m'n  +  1  ij 

and  representing  g(x)  by  a  B-spHne  of  degree  m+n  +  1,  Eq.  (4,  3)  can 
be  written  in  the  form 

£6k'5min  +  l,’‘-k'’t’  =  '4' 51 

k  i  .1 

where  Ax  is  the  sampling  interval.  Equations  (4.3)  and  (4.5)  show 
that  the  B-spline,  which  is  interpo’^.ing  the  deterministic  part  of 
the  blurred  image,  must  be  of  higher  degree  than  the  B-splines 


interpolating  the  object  and  point-spread  function.  In  other  words. 


since  the  blurred  image  is  always  smoother  than  the  object,  a  higher 
degree  spline  can  follow  the  image  function  better  than  one  approxi¬ 
mating  the  object  function.  This  can  be  explained  in  the  Fourier 

th 

domain  by  observing  that  the  Fourier  transform  of  a  m  degree 
B-spline  is  a  sine  function  to  the  power  m+1.  As  m  increases  the 
amplitude  of  higher  frequencies  decreases.  Since  a  blurred  image 
has  less  higher  frequency  content  than  the  object,  a  higher  order 
B-s  pline  can  represent  the  image  better  than  the  one  representing 
the  object. 

Using  vector  space  notation,  Eq.  (4.5)  may  be  written  as 

£  =  Hi  (4.6) 

where  £  and  f_are  vectors  consisting  of  coefficients  g^  and  f  ,  and  hi 

is  a  circulant  matrix  with  elements  h  .  If  the  point-spread-function 

is  of  finite  width,  the  matrix  H  is  banded. 

As  an  experiment  to  compare  spline  functions  with  the  pulse 

approximation  method,  a  rectangular  object  is  blurred  analytically  by 
th 

a  4  degree  polynomial  of  the  form 

“  )  -5- 5S*S5- 5  >4-7' 

=  0  ,  elsewhere 

and  this  is  plotted  in  Fig.  4-1.  The  object  is  a  rectangular  function, 
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ORIGINAL — H  /- - BLURRED 


therefore  it  is  i nti' rpolatod  by  B-splines  of  degree  7. ero. 


!he  so  ond 


derivative  of  h  at  points  x  -  -3.  5  and  x  =  3.  5  is  a  step  fum  tion  and  it 
is  interpolated  by  H -splines  of  degree  two.  Since  the  convolution  of  a 
zero  degree  and  a  second  degree  B-spline  is  a  cubic  spline,  the 
image  is  interpolated  by  cubic  B-splines.  Figure  4-lb,  the  restored 
image  with  and  without  splines,  shows  that  the  spline  restores  the 
edges  much  sharper  and  generates  less  undulations  than  the  common 
pulse  approximation  method.  Using  different  degrees  of  B-splines 
for  object,  image  and  point-spread  functions  depending  on  their 
characteristic  has  led  to  a  more  accurate  model  and  thus  a  better 
quality  restoration.  Spline  restoration  can  also  be  applied  to  a  two- 
dimensional  blur  with  very  good  results  [43]  . 

4.  2  Restoration  of  Space  -  Variant  Degradations  by  the  Minimum 
Norm  Principle 

In  the  previous  chapter,  a  noiseless  blurred  image  was  modeled 
by  the  expression 

£  =  H£  (4.8) 

where  _H  represents  the  blurred  matrix.  If  H  is  square,  non- 

A 

singular  and  well-conditioned,  the  restored  image  f  can  be  obtained  by 

f_  =  H_1g.  (4.9) 

where  H  denotes  the  inverse  of  H.  In  most  practical  cases,  H  is 

either  singular  or  ill-conditioned  due  to  its  large  size  and  due  to  the 

fact  that  most  imaging  systems  irreversibly  remove  certain  aspects 
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of  the  original  object. 


In  underdetermined  or  ovc  rdeter  ri'.i  tied  models. 


which  will  be  defined  later,  H  is  not  a  square  matrix.  Thus  even  in 
the  absence  of  noise  an  estimate  of  f_  cannot  be  obtained  by  '4,  (U. 
This  suggests  the  definition  of  a  reasonable  fidelity  criteria  which 

A 

leads  to  a  unique  solution  for  _f.  The  minimum  norm  criterion  is 
defined  as  the  following 

2 

minimize  )'  f  l1  (4.10) 

among  all  feR  which  minimizes 

!I&-Hfj|2  (4.11) 

where  | 1  •  ! !  denotes  L^-norm  of  the  vector.  Albert  [44]  has  shown 

that  there  exists  a  unique  solution  for  the  above  minimization 

problem.  The  solution  to  (4.  10)  and  (4.  11)  may  be  obtained  by  the 

standard  methods  of  the  calculus  of  variations.  Using  Lagrangian 

2 

parameter  6  ,  the  functional 

W(f_)  =  II*- H  fj!2  +  62!1  fjl 

must  be  minimized.  Taking  a  derivative  with  respect  to  f_, 

=  -2Ht(g-H  f )  +  262f  =  0 
of  —  — 

the  optimal  estimate  for  f_  is 

£_  =  lim  (H^H  +  S2!)"1!^*  (4.12) 

6-*0 

where  I  represents  the  identity  matrix  whose  dimensionality  is 
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understood  from  the  context.  For  example,  if  H  is  m  X  n,  I  is  an 
n  x  n  identity  matrix.  It  is  shown  [44]  that  for  any  m  x  n  matrix  H 

H+  =  lim  t  62lf ’h1  (4.13) 

6-*  0 

always  exists  and  H.+  is  called  the  pseudoinverse  of  H.  For  any  m 
dimensional  vector 

f_  =  H+£  (4.14) 

is  the  vector  of  minimum  norm  among  those  whi 'h  satisfy  (4.  11). 

*  t  t 

The  minimum  norm  f_  is  an  element  of  ^(14  ),  the  range  of  H  ,  and 

satisfies  the  relation 


A 

where  g  is  the  projection  of  £  on  ^>(H).  Since 

(HW  +  62Ht)  =  Hf(H  Hl  4  62I)  =  + 

t  2  t  2  2 

and  since  (H_  H  +  6  I)  and  (H  H  4  6  1)  have  inverses  when  6  >0,  it  is 

clear  that 

(hSi  +  62D~1Ht  =  Hfc(H  Hfc  4  62I)_1 
4* 

and  14  can  also  be  expressed  as 

H+  =  lim  Hfc(H  Hl  4  6ZI)_1  (4.15) 

6-*0 

Equivalently,  the  pseudoinverse  of  a  m  x  n  matrix  14  is 
defined  as  an  n  x  m  matrix  X  satisfying  the  following  four  properties: 
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(4. 16a) 


1) 

H  X  H 

=  ii  . 

2) 

X  H  X 

=  X. 

3) 

,'H  X  )fc 

=  HX 

4) 

(X  H)4 

=  X  H 

(4.  16b) 
(4.  1 6c) 
(4.  1 6d) 


The  above  properties  are  necessary  and  sufficient  conditions  for 
X  =  H*  given  by  (4.  13)  or  (4.  15). 

When  object  and  image  are  represented  by  other  basis  functions, 
such  as  B -splines,  in  a  continuous -continuous  model,  a  similar 
minimization  criterion  may  be  applied.  Let 


f(?) 


g(x) 


m  i 


t 


giBn(x-xi) 


(4. 17a) 


(4. 17b) 


where  and  B^  are  B-splines  of  degrees  m  and  n.  Here  g(x)  is 

related  to  f(x)  by  the  superposition  integral  given  by  (3,5),  Defining 
the  following  objective  function 

<s  ®  2 

w (£ )  =  J  (g<*>  -j  h{x,  §)f(P)d^  dx  +  62j  (f(?>)2d r,  (4.  18) 

_  00  .9 

substituting  (4.  I7a,b)  into  (4.  18)  and  taking  derivatives  with  respect 
to  {_,  the  optimal  estimate  for  f_  is 

*  2  - 1 

f_  =  (P  +  6  BJ  Q  £  (4.19) 

A 

where  the  vectors  f_,  £,  and  the  matrices  P,  B_  and  (3  are  defined  by 
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_  A  A  +>  f- 

t  =  . fj 


m J 


B 


B 


g.  =  [g^  %2’  *  ’ '  ’  gN^ 

(?)  =  [B  (?-?.),  B  fP-?-) _ ,  B  (?-*  )f 

m  1  m  2  m  M 

(x)  =  [B  ix-x  ),B  (x-x  ) . B  (x-x  )]1 

n  1  n  £  n  N 

00 

p(x)  =  j  h(x,  ?)B  (?)dp 

j  m 

00 

P  =  J  £(x)£t(x)dx 

.00 

03 

B  =  f  B  (P)Bt  (P)dP 
—  J  — m  — m 
.00 

CO 

Q  =  J  £(x)JB^(x)dx 

„0D 

The  matrix  _P  is  symmetric  and  non-negative  definite.  Assuming  q 
to  be  an  arbitrary  vector  of  dimension  M,  then 

00  oo 

t  r  t  t  r  t  2 

q  pq  =  J  q  £(x)£  (x)qdx  =  j  (q  £(x))  dx  2  0. 


The  matrix  B  is  symmetric,  banded  and  positive  definite  matrix 

consisting  of  the  values  of  a  B -spline  of  degree  2n-4 1  at  its  knots. 

2 

Therefore,  (P+6  B)  is  positive  definite  and  invertible. 

If  the  functions  f_,  £  and  h  are  represented  by  B -splines  of  degree 
zero,  then  P^  =  H^H,  =  H 1  and  B_  =  I,  and  FTq,  (4.  1  9 >  simplifies  to 
(4,12).  A  similar  formula  is  derived  for  the  continuous-discrete 
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model  when  the  objective  function  is  defined  by  the  minimization  of 


the  second  derivative  of  f  [45]  . 


The  specific  structure  and  properties  of  a  matrix  are  quite 
useful  in  determining  its  pseudoinverse.  For  a  nonsingular  square 

A  _  J 

matrix,  since  f  =  Id  j^is  the  only  vector  minimizing  (4.  11),  the 
pseudoinverse  is  the  same  as  the  inverse.  If  matrix  Id  is  diagonal: 
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where  H+  is  defined  by  (4.21). 

Equation  (4.22)  shows  the  radically  discontinuous  nature  of 
pseudoinversion.  Two  matrices  may  be  very  close  to  each  other 
element  by  element,  but  their  pseudoinverses  differ  greatly.  For 
example,  the  diagonal  matrices 


are  close  to  each  other,  but 


differ  greatly.  The  reason  is  that  (4.22)  exhibits  an  infinite  dis¬ 
continuity  at  X  =  0.  This  characteristic  induces  serious  computa¬ 
tional  difficulties,  particularly  due  to  computer  precision  and  round¬ 
off  error  where  a  small  number  might  be  actually  zero  or  vice  versa. 
This  will  be  discussed  more  in  the  computation  of  H  for  a  general 
m  X  n  matrix. 

The  pseudoinverse  of  a  symmetric  matrix  can  be  derived  by 
using  the  diagonalization  theorem.  A  symmetric  matrix  H  can  be 
written  as 

H  --  Efl  El  (4.2  3) 

where  E  is  an  orthogonal  matrix  and  A_  is  diagonal.  Substituting 
(4. 23)  in  (4. 12)  gives 
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+  2  2  -1  t 

H_  =  lim  EjA_  +6  I)  AE_ 

6-*0 

Z  2-1  t 
=  E_[  lim  (A_  +  6  I)  A]  E_ 

6-»0 

=  EA^E*  (4.24) 

where  A_+  is  defined  by  (4.21)  and  (4.22).  Thus,  the  pseudoinverse 
for  a  symmetric  matrix  is  obtained  by  pseudoinver ting  the  diagonal 
matrix  of  its  eigenvalues.  Equation  (4.24)  can  equally  be  expressed  as 

H+  =  X.*eeef  (4.25) 

where  e  is  the  eigenvector  of  H  associated  with  the  eigenvalue  X.. 

— i  —  i 

If  H  is  a  rectangular  matrix  of  full  row  rank,  i.e.,  the  rows  of 
II  are  linearly  independent,  then  H  _h' t  is  invertible  and  eq.  (4.  15) 
simplifies  to 

H+  =  Hl(H  hS"1  (4.26) 

Although  Eq.  (4.25)  presents  a  straightforward  method  for  computing 
H+,  the  problem  of  inverting  the  m  xm  matrix  (H^  Ed*" )  remains.  This 
can  cause  difficulties  for  even  moderate  size  images.  In  this 
situation,  the  observation  can  be  partitioned  into  smaller  segments 
which  are  used  for  estimation  of  the  corresponding  object  sections. 
Moreover,  since  the  number  of  linear  equations  is  less  than  the 

A 

number  of  unknowns  in  (4,8),  the  estimated  object  f_is  not  necessarily 
equal  to  the  original  object  In  other  words,  a  full  recovery  of  the 
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vector  i_  is  not  guaranteed.  Nevertheless,  the  estimated  vector  {_ 
satisfies  the  relation 

H  f_  =  £_  =  H  L 

because 

H  f_  =  H  H+g.  =  H  Hl  ( H  H  )  £.  =  £.• 

If  H  is  a  rectangular  matrix  of  full  column  rank,  i.e.,  the 
columns  of  H  are  linearly  independent,  then  ^His  invertible  and 
Eq.  (4.  15)  simplifies  to 

H+  =  (HWH  <4*27> 

In  this  case,  the  number  of  linear  equations,  or  in  other  words, 
the  number  of  observations,  is  more  than  the  number  of  unknowns. 

If  the  system  of  equations  is  consistent,  i.e.,  £  is  a  linear  combin¬ 
ation  of  column  vectors  of  H,  which  is  always  true  in  noiseless 
models,  the  estimated  object  £_  is  the  same  as  the  original  object  f 

because 

1  =  H+£  =  (hWV*  =  (HWVHf  =  i_  . 

Thus  the  object  can  be  recovered  from  the  image  without  any  error. 
A  computationally  efficient  method  utilizing  the  Fourier  properties 
of  circulant  matrices  has  been  introduced  for  the  pseudoinversion  of 
full  column  matrices  in  space -invariant  degradation  [n]  . 

When  H  is  a  general  m  X  n  matrix  with  no  particular  structure, 
singular-value-decomposition  (SVD)  techniques  can  be  used  for 


In  practical  cases,  a  judicious  choice  of  eigenvalue  cutoff  X^ 
must  be  made  for  nonzero  eigenvalues.  If  the  X.'s,  ordered  in 
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decreasing  value,  show  a  sudden  decrease  in  value  as  a  func  tion  of 


the  index  i,  then  the  threshold  may  be  located  at  that  point.  The 
decrease  in  value  can  be  by  a  factor  as  small  as  the  machine 
precision.  If  such  a  sudden  decrease  does  not  exist,  a  threshold  0 
which  is  dependent  on  machine  precision  must  be  selected  and  the 
eigenvalues  smaller  than  e  are  declared  zero.  The  value  r  deter¬ 
mines  the  rank  of  H  and  small  eigenvalues  X  X  are  assumed 

to  be  roundoff  error.  Equation  (4.29)  expands  H  in  terms  of  system 
eigenvectors;  thus  the  X 's  are  the  effective  spectral  components. 
General  outer-product  expansions  of  H  are  given  by 

H  =  ^  0,  .u.  vt  (4.31) 

1  J 

where  u.  and  v.  can  be  the  discrete  Fourier  basis  vectors,  Walsh- 
1  J 

Hadamard,  Haar,  Slant,  or  other  orthonormal  bases.  With  a 

space-invariant  degradation,  H  is  a  circulant  matrix  that  can  be 

diagonalized  by  discrete  Fourier  transforms.  Thus,  the  SVD 

procedure  is  analogous  to  the  discrete  Fourier-inverse-filtering 

method  that  is  widely  used  for  space -invariant  processing. 

4.4  Restoration  of  Astigmatism  and  Curvature-of-Field 

Optical  images  are  subject  to  a  number  of  blurring  effects  due  to 

aberrations.  Certain  aberrations,  such  as  spherical  aberration, 

can  be  described  by  convolution  integrals  and  canbe  solved  in  the 

Fourier  domain.  For  other  aberrations  such  as  coma,  astigmatism 

55 


or  curvature-of-ficld,  the  blurring  is  spare-variant.  When  the 
effects  of  astigmatism  and  cur vature -of-fi eld  predominate,  the 
geometrical-optics  aberration  functions 


x-  “ 


(2C+D)-  r  cos 


(4.  32a ) 


y-P  = 


r  sin  0 


(4.  32b ) 


describe  the  displac  ement  of  an  image  point  from  its  ideal  (Gaussian) 
intercept  in  the  image  plane.  Here  r  and  6  are  ray  intercepts  in  the 
exit  pupil  of  the  optical  system,  and  C  and  D  are  constant  coefficients 
describing  the  degree  of  astigmatism  and  curvature -of-field, 
respectively  [.22]  ,  [48*1  .  Using  a  technique  described  in  [22]  and  [48], 
the  space -variant  point-spread  function  (SVPSF)  of  the  system  for  a 
circular  exit  pupil  of  radius  R  is  obtained  as 


h(x,  y;ff,  T]=0)  = 


1 


(x-») 


D(2C  +  D): 

0  . 


2  2  4 

DR? 


elsewhere 


2  2  4 
(2C -t-D)  R  c 


<;  1 


(4.  33) 


assuming  an  object  impulse  function  at  (?,  TbO).  This  function  is 

given  in  Fig.  4-2  for  the  impulses  at  various  locations  in  the  ( F ,  T) ) 

plane.  The  region  of  nonzero  response  are  defined  by  ellipses  which 

increase  in  size  proportional  to  the  square  of  the  radial  distance,  and 

4 

the  amplitude  of  the  response  decreases  inversely  with  c  .  Although 
the  system  is  strongly  space-variant  and  the  blurring  occurs  in  both 
radial  and  angular  directions,  changes  in  the  amplitude  and  shape 
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of  the  response  are  a  function  only  of  radial  distance.  Thus,  the 
PSF  h(*  )  of  Eq.  (4.  33)  written  with  7)  =  0  is  just  rotated  about  the 
origin  to  obtain  the  general  response.  Because  of  the  inherent 
circular  symmetry,  the  system  complexity  can  be  reduced  by  a  polar 
coordinate  transformation  of  the  form 

f  =  p  cos  CPq  (4.  34a) 

7)  =  Oq  sin  ©q  (4.  34b) 

in  both  object  (with  subscript  0)  and  image  (without  subscript) 
coordinates,  and  rewriting  Eq.  (4.32)  in  the  form 

c  -  Oq  =  (2C  +  D)0gr  cos  9  (4.  35a) 

Cp  -  CPq  =  tan  [Dp^r  sin  9/  [l  +  (2C  +  D)Dor  cos  91}  =  u(p  ) 

(4. 35b) 

where  (p  ,  CO  )  and  fn.cp)  are  the  object  and  image  polar  coordinate 
variables.  In  this  form,  the  two-dimensional  space-variant  radial 
blur  becomes  decoupled  from  the  angular  blur  because  Eq.  (4.  35a) 
does  not  contain  cp  or  rp^.  The  blur  in  the  angular  direction  is  space- 
invariant  in  CD  and  a  slowly  varying  function  of  position  as 
expressed  by  ufo^)  in  Eq.  (4,  35b). 

When  the  degradation  is  purely  astigmatic  with  no  curvature -of- 
field,  the  D  coefficient  in  Eq.  (4.  32)  becomes  zero  and  Eq.  (4.  33) 


becomes  singular  because  no  blurring  occurs  in  the  angular  direction. 
To  find  the  SVPSF  for  astigmatism  only,  Eq.  (4.  33)  is  first  collapsed 


r 


| 


to  a  purely  radial  space -variant  blur  by  h^fx;^,  T)  =  0)  by  evaluating 

00 

ha(x;F,Ti=0)  =  J  h(x,  y;F,T)=0)dy  (4.36) 

.00 

and  taking  the  limit  as  D  approaches  zero.  The  result  is 
r  2  2  4  2-,y 

h  (x;*, T)=0)  =  - ^ - L- J—  ,  ?-2CRf  Sx<Ft2CR?  (4.37) 

a  '  2cV 

and  zero  elsewhere.  The  region  of  support  and  cross  section  of  this 
function  are  shown  in  Fig.  4-3.  With  astigmatism  only,  the  degrada¬ 
tion  reduces  to  a  two-dimensional  space-variant  line  blur  in  a  purely 
radial  direction.  Figure  4-4a  is  an  aerial  photograph  displayed  as 

128  a  128  discrete  picture  elements  after  blurring  by  astigmatism 

-4 

with  R  =  1  and  C  =  7.  5  x  10  .  Note  that  the  blurring  increases  from 

zero  on  the  optical  axis  (upper  left  corner)  to  nearly  50  picture 
elements  in  width  as  a  function  of  increasing  radius. 

For  the  system  degradations  due  entirely  to  astigmatism,  the 
ideas  of  coordinate  transformation  restoration  (CTR)  [22"]  can  be  used 
with  little  modification.  The  basic  idea  is  to  reduce  space-variant  to 
space -invariant  distortions  by  invertible  coordinate  transformation.  The 
SVPSF  of  Eq.  (4.  37)  can  be  modeled  as  a  polar  coordinate  transform¬ 


ation  on  the  object  coordinates,  followed  by  identical  space-variant 
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(a)  SVPSF  for  pure  astigmatism  with  inputs  at  various 
distances  f , 


Figure  4-3.  SVPSF  and  its  cross  section  for  pure  astigmati 


s  m. 
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fa)  Image  degraded  by 
astigmatism 


(b)  Polar  transformation  of  (a) 


Figure  4-4.  Restoration  of  astigmatism. 


radial  blurring  for  each  cp  variable.  A  final  inverse  polar  trans¬ 
formation  on  the  image  coordinates  completes  the  model  and  produces 
the  PSF  shown  in  Fig.  4-3.  Performing  this  decomposition  reduces  a 
four -dimensional  space-variant  restoration  problem  to  a  single  two- 
dimensional  problem. 

Using  the  transformation,  we  can  write  the  radial  degradation 
in  matrix  form  as 

G(p,tp)  =  H(d,d0)F(o0,V)  (4.381 

where  and  F^  are  matrices  representing  image  and  obioct  in  polar 
coordinates  and  Id  is  the  blur  matrix  obtained  by  applying  a  mono¬ 
spline  quadrature  formulae,  which  was  discussed  in  section  3.5,  to 
the  continuous  space  description  ofEq.  (4,371.  This  quadrature 
formulae  provides  a  more  accurate  and  smooth  discrete  approxima¬ 
tion  to  the  continuous  representation  of  Eq.  (4.  371. 

The  space-variant  restoration  procedure  (CTR)  proceeds  by 
inverting  the  two  polar  coordinate  distortions  and  solving  Fq.  (3.  381. 
Unfortunately,  a  direct  inversion  of  H_  is  usually  not  possible  because 
point-spread  function  matrix  H  tends  to  be  ill-conditioned  leading  to 
numerical  problems.  The  ill-conditioning  is  a  result  of  the  inform¬ 
ation  loss  associated  with  the  imaging  process;  thus  H  is  generally 
singular  and  pseudo-inversion  must  be  used.  For  inversion  of 
Fq.  (4.38),  the  singular  value  decomposition  algorithm,  which  was 
discussed  in  the  previous  section,  is  used  to  obtain  a  unique  pseudo- 


inverse  which  is  then  used  in  the  restoration  operation 

F(p0,cp)  =  H  +  (P0,p)G(p,C?)  (4.39) 

The  C  TR  procedure  has  been  implemented  or.  the  imag'  degraded 
only  by  astigmatism  in  Fig.  4-4a.  First  a  polar  coordinate  trans¬ 
formation  is  performed  to  produce  Fig.  4-4b  in  which  the  space- 
variant  blur  (4.  38)  occurs  in  only  the  radial  direction.  Figure  4-4c 

shows  the  restoration  by  SVD,  in  which  the  7  singular  values  out  of 

_  5 

128  whose  magnitudes  were  less  than  10  were  not  used.  Figure  4-5 
shows  the  singular  values  in  a  decreasing  order.  Following  restor¬ 
ation  of  each  line  in  the  (p,Cp)  system,  an  inverse  polar  coordinate 
transformation  is  used  to  produce  the  final  result  of  Fig.  4-4d. 

This  procedure  can  also  be  used  for  restoration  with  both 
astigmatism  and  curvature-of-field  present.  First,  the  imaging 
equation  is  expressed  with  a  polar  coordinate  transformation  (3.34) 
in  the  form 

CD 

g(p,cp)  =  II  h(P-5P:Po'5P0)f(Po’'co)dpod(:po  (4.40) 

-m 

and  then  rewritten  as 

00 

g(p»<P)  =  JJ  h'(p,  D0,cp-cp0,  u(p0))f(p0>cp0)dp0dcp0  (4.41) 

_  CD 

using  CP-CPq  and  u(o^)  of  Efcj.  (4.  35)  to  emphasize  the  functional 
dependence.  Defining  a  Fourier  transform  of  g(nfp)  in  the  ©variable  by 
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(4.42) 


r 


g(p,X)  =  J  g(P,  cp)exp(-j2nXcp)dcp 

>00 

the  transform  of  both  sides  of  Eq.  (4.41)  is  taken  to  obtain 

CO 

g(p,X)  =  jj  f(o0,«PQ)H(D,  P0»  X.  u(o0))exp(-j2rrXcp0)dcp0dp0 

(4.43) 

where  h  is  the  Fourier  transform  of  h'  in  cp.  Grouping  terms  con¬ 
taining  cp  on  the  right  side  of  Eq.  (4.43)  enables  a  transform  in  this 
variable  to  be  evaluated.  The  resulting  transformed  function  f  (p,  X) 
is  given  by 

so 

f(o0.M  =  J  f(o0,CP0)exp(-j2TTXcp0)dcD0  (4.44) 

>QD 

and  the  reduced  system  equation  obtained  from  Eq.  (4,41)  is 

oo 

g(D,X)=J  h(P,pQ,  X)f  (pQ,  X)dpQ  (4.45) 

.00 

where 

Mo.Pq.X)  =  h(p,p  ,  X,u(pQ)) 

is  written  as  a  function  of  three  variables  to  show  explicit  depend¬ 
ence. 

This  procedure  can  be  extended  for  the  restoration  of  images 
degraded  by  simultaneous  astigmatism  and  curvature -of-field 
aberrations.  Following  a  polar  coordinate  transformation,  a  Fourier 
transform  in  cp  as  expressed  by  Eq.  (4,42)  is  performed  to  partially 
decouple  a  blur  as  a  slowly  varying  function  of  u(Pq).  The  reduced  ^ 


System  given  by  the  continuous  space  integral  Eq.  (4.45i  has  the 
same  form  as  the  discrete  equation  (4.  X8)  for  astigmatism.  An 
estimate  ffn^.X)  is  then  produced  by  the  SVD  for  each  separate  X 
using  similar  techniques.  The  computational  effort  in  this  operation 
may  be  reduced  by  using  the  known  variation  of  Mo.Pq.X)  with  X  . 

After  the  entire  f(p^,X)  has  been  obtained  a  series  of  one -dimensional 
inverse  Fourier  transforms  in  cp^  is  taken  to  find  f ( ,  Cp^ ) ,  and  an 
inverse  polar  coordinate  distortion  is  used  to  get  f(~,T)  as  the  final 
restored  object.  This  procedure,  while  requiring  large  capabilities 
in  computing  and  storage,  is  the  only  practical  method  for  restoration 
of  images  of  even  moderate  size.  The  general  four -dimensional 
space-variant  blur  is  effectively  reduced  to  a  set  of  space-variant 
two-dimensional  problems  whose  point-spread  functions  depend  in  a 
weil-known  way  on  X  . 

4.  5  Overdetermined  and  Underdetermined  Models 

When  the  degradation  matrix  Id  of  an  imaging  system  is  of  full 
column  rank,  the  model  is  called  overdetermined.  In  practice,  this 
usually  occurs  in  two  situations.  The  first  is  when  the  object  has 
zero  background  and  the  object  and  image  are  sampled  at  the  same 
rate.  The  second  occurs  when  the  image  is  sampled  at  a  higher  rate 
than  the  object. 

Suppose  the  object  function  f  has  zero  background,  i.e.,  f(x)  is 
zero  if  x  <  Xj  ir  x  >x^,  and  the  point-spread  function  h(x)  is  space- 
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invariant,  symmetric  and  space -limited  of  width  2.1,  Then  the  limits 
for  the  convolutional  integral  are 


a  =  max(x-Jt,  Xj)  (4.46a) 

b  =  min(x+ 1,  x^ )  (4.46b) 

and  the  degraded  image  is  given  by 

b 

gfx)  =  J  h(x-ff)f(f )dP  .  (4.47) 

a 

Assuming  uniform  sampling  of  image,  object  and  PSF  with  sampling 
interval  Ax  =  1,  the  continuous  model  can  be  discretized  as 


g(t) 


c.  ,h(i-j)f(j) 


(4.48) 


where  c.  .  is  the  quadrature  coefficient  associated  with  h(i-i),  and 
'-J 

K  j  =  max(t-L,  1 )  (4.49a) 

K2  =  min(i+L,  N)  (4.49b) 

where  L,  is  the  integer  part  of  l.  The  image  sample  g(i)  is  not  zero 
if  -L  +  l  ^N+L,  and  thus  the  number  of  observations  is 

M  =  N  +  L  -  (-L+I  )  +  I  =  N  +  2L  (4.50) 


Using  a  vector  space  notation,  the  equation  (4.48)  becomes 


£  =  Hf 


(4.  51) 
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where  ff  is  the  overdetermined  M  X  N  blur  matrix  defined  by 


since  f_  j  -_f  ^  ^  0>  Eq.  (4.54)  indicates  that  a  linear  combination  of 
columns  of  H  is  zero  which  contradicts  the  assumption  of  an  over¬ 
determined  model.  This  unique  solution  can  be  obtained  by  the 
pseudoinverse  of  a  full  rank  matrix  of  by  SVD  algorithm  discussed 
in  section  4. 3.  In  noisy  images,  ^  is  not  necessarily  in  R  (HJ,  and  a 
least  square  criteria  based  on  (4.  11)  leads  to  a  unique  solution  [161  . 

A  more  realistic  model  for  an  imaging  system  can  be  obtained  if 
no  restrictions  are  imposed  on  the  background  of  the  object.  In 
practice,  few  objects  are  recorded  with  a  background  of  zero  or 
known  intensity.  Moreover,  because  of  computational  problems,  an 
image  is  often  partitioned  into  sections  before  being  processed,  and 
thus  the  assumption  of  zero  background  cannot  be  valid.  A  model 
must  be  used  that  relates  a  portion  of  the  object  to  the  corresponding 
segment  of  the  image  without  any  restrictions  on  the  background.  In 
such  a  model,  because  of  blurring,  a  portion  of  the  image  is  affected 
by  a  larger  segment  of  the  object.  Therefore,  if  the  object  and  image 
are  sampled  with  the  same  rate,  the  matrix  _H  has  more  columns 
than  the  rows,  and  the  system  is  underdetermined.  Following  the 
same  procedure  as  with  overdetermined  models,  the  system  is 
represented  by 

£  =  H£  (4.55) 

where  _H  is  M  x  N  blur  matrix,  and 

(><> 


M  =  N  -  2L 


The  matrix  H  is  given  by 

H  = 


where  h  and  c_  are  the  same  as  overdetermined  model.  As  mentioned 
in  section  4.  3,  the  system  (4.55)  does  not  have  a  unique  solution. 

The  minimum  norm  solution  based  on  (4.  10)  and  (4.  11)  is  unique  and 
can  be  obtained  by  (4.26)  if  H  is  of  full  row  rank,  or  by  the  SVD 
algorithm. 

4.  6  Experimental  Results 

To  illustrate  image  restoration  by  pseudo-inversion,  Fig.  4-6a  is 
selected  as  a  test  scene  which  is  originally  of  size  128  x  128  picture 
elements  'pixels)  with  1  1  0  X  1 1  0  nonzero  elements.  For  di splay ,  this 
zero  background  picture  has  been  enlarged  by  cubic  spline  interpolation 
to  an  image  of  size  256  x  256  pixels.  Figure  4-6b  represents  the 
image  after  undergoing  motion  degradation  with  a  blurring  point 
spread  function  of  length  11  pixels.  This  is  a  severe  blur  because 
it  is  1  /  1  0  of  the  original  picture  size.  J'he  observed  image  is  of 
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(4.  56) 
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(b)  Blurred  Icl  Restored 


Figure  4  b.  Restoration  of  motion  degradation,  overdetermined 
model. 
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size  1  10  /  120  and  ihus  t  lu*  system  is  Dvi-rdcarmincd  as  disv-'isscd  in 
sec  -  i or.  4.5.  I' he  blur  mal  r  i  x  H_i  s  a  1  2 0  x  1  1  0  ba  nd ed  ma  t  r  i  x  g  i  ve  n 
by  Fq.  (4.  52)  with  1  1  nonzero  elemen  s  on  each  column.  All  hough 
this  matrix  is  of  full  column  rank,  because  of  high  dimensionality 
and  ill -conditioning,  using  Fq.  (4.27)  and  inver.ing  hSi  produces 
computational  difficul  ies.  Thus,  the  SVD  algorithm  has  been 
employed  to  compute  the  pseudo-inverse  of  matrix  H.  Figure  4-7 
shows  the  singular  values  of  H  ordered  in  decreasing  value.  Since 
the  singular  values  are  much  greater  than  the  machine  precision  and 
there  is  not  a  sudden  decrease  in  value,  all  of  them  have  been  used 
used  for  computation  of  H+.  Figure  4-6c  shows  the  restored  image. 
As  far  as  visual  perception  is  concerned,  this  restored  image  is 


identical  to  the  original  scene. 


Chapter  5 

RESTORATION  OF  NOISY  IMAGES 

In  the  previous  chapter,  the  problem  of  image  degradation  and 
restoration  was  modeled  for  a  noiseless  system.  In  reali'y,  an 
image  is  often  affected  by  a  variety  of  noisy  sources.  The  scanner 
measurement  error  is  a  source  of  noise  which  adds  some  element 
of  uncertainty  to  the  measured  signal.  The  quantizer,  that  maps  the 
signal  amplitude  to  a  finite  number  of  digital  levels  for  computational 
and  coding  purposes,  is  ano.her  source  of  noise.  Coding  and  channel 
errors  occur  when  he  image  is  transmitted  through  a  noisy  channel. 
Computers  produce  truncation  and  roundoff  errors  due  to  the  limits 
of  machine  precision.  Photographic  film  is  the  most  common 
recording  system  in  image  processing  and  is  another  source  of  noise 
which  will  be  discussed  in  more  detail  in  section  5.  3. 

Although,  these  examples  are  not  an  exhaustive  list  of  noise 
sources,  noise  of  various  types  is  a  common  problem  in  every  type 
of  imaging  system.  In  this  chapter  the  effects  of  noise  in  image 
restoration  is  considered  from  several  viewpoints. 

5.  1  Discrete  Wiener  Filter 

The  basic  idea  of  the  Wiener  filter  and  its  derivation  in  a 
continuous  model  was  discussed  in  section  2.  3.  In  this  section, 
the  same  fidelity  criteria  is  applied  to  a  discrete  model  and  the 
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filter  is  obtained. 
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singular.  If  the  noise  is  white,  this  condition  is  satisfied,  because 


H  C  „  H  is  non-negative- definite  and  C  is  posit  ive -definite ,  there 
- f —  — n 

fore  Id  C^H1  +  C_  is  positive-definite  and  consequently  is  nonsingular. 
Another  version  of  (5.  5)  can  be  obtained  by  using  the  maLrix  identity 
[51] 


t  t  -1  -1  t  -1  -1  t  -1 

A  B  (C  +  B  A  B  )  =  (A  +BC  B)  B  C 


(5.6) 


in  B}.  (5.  5) 

t  - 1  - 1  t  - 1 

R  =  (H  C  H  +  C,  )H  C  (5. 7) 

— - n  —  — f  —  n 

where  C  and  C  are  assumed  to  be  nonsingular.  Although  these  two 
— n  — i 

versions  are  equivalent,  their  computation  in  practice  depends  on  the 

structure  of  H,  C  and  C,.  If  H  has  fewer  columns  than  rows, 

—  — n  — f  — 

H  C-H^  P  C  has  a  higher  dimension  than  F^C  +  C,  ,  and  there  - 
fore  needs  more  compulations  for  inversion.  The  order  is  reversed 
when  the  number  of  columns  is  more  than  the  number  of  rows. 
Moreover,  Eq.  (5.7)  requires  the  inversion  of  two  more  matrices; 

C  and  C, . 

The  Wiener  filter  can  also  be  obtained  by  minimization  of  the 
objective  function 

W(n  =  II  C^ill2  +  ||Cn?(g.-H_f)||2  (5.8) 

l  X 

where  C ^  2  and  C^2  are  the  whitening  filters  for  the  signal  and  noise, 
respectively.  Faking  derivatives  with  respect  to  f_ 
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dW 

df 


C^f  +  HfcC  1  H  f  -  H^C ' 1  g 
— f  — - n - - n  h- 


and  setting  the  result  equal  to  zero,  the  estimate  is 


-  t  -1  -1  -1  t  -1 

f  =  (H  C  H  +  C,  )  HC  g  (5.9) 

—  - n  —  — f - n 

which  is  the  same  as  (5.7).  The  inversion  of  large  size  matrices  in 
Eqs.  (5.  5)  and  (5.  7)  is  not  an  easy  task,  particularly  because  of  ill- 
conditioning.  Using  the  Fourier  transform  for  the  inversion  of 
circulant  matrices  leads  to  a  fast  Wiener  filter  with  much  less 
computation  and  higher  efficiency  [13]  . 

5.  2  Filtering  of  Unblurred  Noisy  Images  by  Smoothing  Spline 

Functions 

As  discussed  in  Chapter  1,  the  Wiener  filter  has  many  limitations 
and  shortcomings.  Not  only  does  it  require  the  most  a  priori  inform- 
ation,  it  often  produces  an  estimate  of  the  restored  image  with  poor 
visual  quality  in  comparison  to  other  filters  [52  J  .  l’he  constrained 
least  squares  estimate  which  has  been  introduced  as  an  alternative 
to  the  Wiener  filter  is  capable  of  producing  images  with  much  higher 
visual  quality.  Moreover,  this  estimate  requires  less  a  pr i ori 
information  than  the  Wiener  filter  [20 ]  . 

In  this  section,  a  constrained  least  squares  estimate  based  on  the 
minimization  of  the  second  derivative  and  local  statistic  s  of  the  noise 
is  defined  and  the  corresponding  filter,  which  produces  a  cubic  spline- 
function  as  the  estimate,  is  obtained.  The  parameters  of  this  filter 
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determine  the  local  smoothing  window  and  overall  extent  of  smoothing, 
and  thus  one  can  control  the  tradeoff  between  resolution  and  smoothing 
in  a  spatially  nonstationary  manner.  Since  the  derivation  of  the  filter 
for  unblurred  noisy  images  is  simpler  and  the  estimate  can  be  compu¬ 
ted  by  highly  efficient  algorithms,  the  derivation  of  the  filter  for 
unblurred  noisy  images  is  covered  here.  The  corresponding  filter  for 
noisy  blurred  images  will  be  covered  in  section  5.  5. 

For  an  unblurred  noisy  image,  the  image  gfx)  is  given  by 

g(x)  =  f(x)  +  n(x)  (5.10) 

in  a  continuous  model.  The  additive  noise  n(x)  is  assumed  to  be  an 
uncorrelated  zero  mean  random  process,  i.e.  , 


where  the  positive  number  6.  locally  controls  the  smoothing  window 
at  point  x.  and  S  controls  the  overall  extent  of  smoothing.  If  tin- 
standard  deviation  of  noise  or  its  estimate  at  point  x  is  available,  then 

it  can  be  used  for  6..  In  this  case,  natural  values  of  S  lie  within  the 

i 


7  K 


confidence  interval  of  the  left  side  of  (5,  13)  that  is 


I  i 

N  -  (2N)2  SN  +  (2N)2  (5.14) 

where  N  is  the  number  of  data  points.  Reinsch  [53]  has  shown  that 

the  solution  to  (5.  12)  and  (5.  13)  is  a  cubic  spline  and  more  generally 

is  a  spline  function  of  degree  2K-1  for  least  squares  minimization 
tVi 

of  the  K  derivative  instead  of  the  second  derivative  [54]  .  The  case 
K  =  2  leads  to  very  simple  algorithms  for  the  construction  of  f(x). 
Moreover,  cubic  splines  give  satisfactory  results  and  are  easy  to 
evaluate.  Choosing  S  equal  to  zero  implies 

f(x.)  =  g(x. )  i  =  1 . N 

which  leads  to  the  problem  of  interpolation  by  cubic  spline  functions. 

Applying  the  well-known  Lagrange  multiplier  method,  along  with 
the  auxiliary  variable  z  to  change  the  inequality  constraint  (5.  13)  to  an 
equality  constraint,  the  objective  function 

XN  2 

(f"(x))  dx  +  p 

*1 

must  be  minimized.  The  optimal  function  ffx)  is  determined  as  the 
following 


£< 


g(x  )-f(x. )  2 


6. 


+  z  -S 


(5.  15) 


fix.)  -  ffx. )  =  0  i  =  1 . N  (5.  16a) 

\  -  i  + 

f'fx.)  -  f'fx. )  =0  i  =  1 . N  (5.  16b) 

l  -  i  4 

f"fx. )  -  f'fx.  )  =  0  i  =  1 . N  (5.  16c) 

1  -  1  i 
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(k) 

where  f'  (x) 


f<x.  )-g(x. ) 

f *"(*.)  -  f'"(x. )  =  2p - l— - 1-  i  =  1 . N 

1  -  IT 

i  (5.  1 6d) 

(k) 

lim  f  (x±h)  for  K  =  0,  . .  . ,  3.  Moreover 
h-+0 


f""(x)  =  0,  x.  <x  <x  i  =  1 . N-l.  (5.  1  6e) 

i  1 1  i  ' 

Equations  (5.16)  indicate  that  the  function  f(x)  is  composed  of  piece- 

wise  polynomials  of  degree  3  or  less  in  each  interval  [x  ,  x  *|  of  the 

i  i  +1 

form 

2  3 

f(x)  =  a.  +  b.  (x-x. )  +  c.  (x-x. )  +  d.  (x-x. )  (5.17) 

1  l  1  l  l  11 

such  that  they  are  continuous  up  to  the  second  derivative  at  their 
joining  points,  and  thus  the  solution  is  a  cubic  spline.  Assuming 
f(x)  to  be  a  natural  spline  [33]  of  degree  3,  the  following  extra 
conditions 

f”(x1)_  =  f"'(x1)_  =  f»(xN)+  =  f'"(x  )  =  0  (5.  18) 

must  be  satisfied.  Substituting  (5.  17)  in  (5.  16)  and  (5.  18),  the 
optimal  filter  with  respect  to  conditions  (5.  12)  and  (5.  13)  is  obtained 


as  [Appendix  A] 

A  2  t  2  -It 

!=  g.  -  D  Q(Q  D  Q  +  p_T)  Q&  C>.  IP) 

where 

t 

1  =  lffx1),f(x2) . f(x  )]  (5.20a) 

g.  =  [g(x|),g(x2) . g(xN>l  (5.20b) 
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(5. 20c) 


D  =  diag[61,fc2 . &N]. 

Matrix  _T  is  a  positive  definite  tridiagonal  matrix  of  order  ]\'-2,  and 
Q  is  a  tridiagonal  matrix  with  N  rows  and  N-2  columns.  For  a 
uniformly  sampled  image  with  unit  sampling  interval,  Q  and  T  have 
the  form 


(5.21a) 


4/3  1/3 


1/3  4/3  1/3 

\  X  N. 


1/3  4/3  1/3 


1/3  4/3 


(5.21b) 


Here  matrix  represents  a  second  order  differentiation  matrix. 

Since  Q  is  a  matrix  of  full  column  rank,  and  ID  is  nonsingular,  the 
t  2 

matrix  Q  D  Q  +  p_T  has  an  inverse  for  all  values  of  p  i  0,  Thus,  if  p 

A 

is  known,  the  estimate  f_  can  be  obtained  by  (5.  19).  The  Lagrangian 


parameter  p,  like  S,  controls  the  overall  extent  of  smoothing  and 


may  he  determined  in  terms  of  S. 

The  objective  function  (5.  15)  has  to  be  minimized  also  With 
respect  to  z  and  p,  leading  to  the  conditions 


pz  =  0 


(5.  22) 


f(x.  )-g(x. ) 
i  l 

6. 


) 


Substituting  (5.19)  in  (5.23)  gives 


F(p) 


D  QfQ^Q  +  pn^Q^H  =  (S 


2 

-z  ) 


(5.23) 


(5.24) 


Condition  (5.22)  implies  either  p  =  0  or  z  =  0.  The  first  case  is  only 
possible  if  F(0)  ^S,  and  thus  the  straight  line  fitted  to  data  points  by 
least  square  principles  satisfies  condition  (5.13)  and  the  cubic  spline 
reduces  to  this  straight  line.  If  F(0)  >S,  then  p  t  0  and  z  =  0,  the 
inequality  constraint  (5.  13)  changes  to  an  equality  constraint  and  Eq. 
(5.  24)  can  be  written  as 

F(p)  =  S2  (5.25) 

Reinsch  [40]  has  shown  that  there  exists  a  unique  solution  for  p 
satisfying  (5.  25)  and  (5.  12).  This  positive  unique  solution  can  be 
determined  by  using  Newton's  method. 

Since  the  parameters  p  and  S  are  interrelated  by  (5.25)  and 
both  have  the  same  effect  on  filtering,  by  selecting  p  instead  of  S, 
one  can  skip  the  iterative  Newton's  method  to  compute  p  and  reduce 
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the  numerical  operations  considerably. 

The  matrix  _T  is  positive  definite  and  invertible,  by  using  matrix 
identity  [Appendix  b!  ,  Ea.  (5.  19)  can  thus  be  written  in  the  more 
concise  form: 

- 1  2  - 1  t  1 

i_  =  (I  +  p  D  Q  T  Q  )  £  (5. 26) 

Although  (5.26)  appears  simpler,  it  needs  two  matrix  inversions  in 

comparison  to  Eq.  (5.  19)  which  requires  the  inverse  of  a  banded 

matrix.  The  Cholesky  decomposition  [55]  R_^R  of  a  positive-definite 
t  2 

band  matrix  Q  E)  +  p_T,  where  R_is  a  lower  diagonal  (triangular) 
matrix,  provides  an  efficient  computational  algorithm  for  performing 
the  filtering  ofEq,  (5.19). 

5.  3  Application  of  Smoothing  Spline  Filter  in  Signal-Dependent  Noisy 
Images 

An  interesting  property  of  the  fidelity  criterion  (5.  13)  is  that  the 
smoothing  window  can  be  locally  controlled  by  determination  of  6.. 

If  the  noise  variance  is  higher  in  some  regions,  6.  can  be  set  larger 
at  that  region.  This  property  enables  the  spline  filter  capable  to 
restore  images  even  with  the  difficult  problems  of  signal-dependent  or 
multiplicative  noise.  A  common  source  of  this  type  of  noise  is 
photographic  film,  and  since  it  is  widely  used  as  a  recording  system 
in  image  processing,  the  restoration  of  images  degraded  by  film- 
grain  noise  is  very  desirable. 

Image  formation  on  a  photographic  film  is  a  complex  optical  and 
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chemical  process.  A  very  accurate  model  of  this  process,  even  if 
possible,  is  too  complicated  to  be  used,  and  a  very  simplified  model 
can  hardly  represent  the  actual  physical  process.  A  reasonable  and 
fairly  accurate  model  for  film  grain  noise  is  given  by  r56]  ,  [58] 

1 

D  =  D  +  kD2n  (5. 27) 

r  s  s 

where  D  is  the  signal  density,  is  recorded  density  on  film  and  n 
is  zero  mean  Gaussian  with  unit  variance  that  is  statistically  independ¬ 
ent  of  the  signal.  A  block  diagram  of  this  model  is  shown  in  Fig. 

5-la. 

Equation  (5.27)  shows  that  when  the  signal  has  a  higher  amplitude, 
noise  has  a  greater  variance.  To  apply  the  above  mentioned  filter, 
the  value  of  6.  and  an  estimate  of  the  signal  are  needed.  Hunt  and 

l 

Cannon  [59]  have  shown  that  images  are  well  described  statistically 
as  a  stationary  variance  about  a  spatially  non-stationary  local 
mean.  Assuming  ergodicity  of  similar  classes  of  images,  the 
spatial  average  may  be  used  as  an  estimate  of  the  ensemble  average. 

A  local  spatial  average  is  thus  used  as  a  non-stationary  estimate  of 
the  signal  mean.  Since  a  linear  estimate  of  6.  in  terms  of  the  observ¬ 
ation  is  used  in  a  nonlinear  manner  to  filter  the  image,  the  overall 
filtering  is  nonlinear.  The  block  diagram  of  this  filter  is  shown  in 
Fig.  5-lb. 

The  filter  of  Eq.  (5.  19)  has  been  applied  to  an  image  corrupted 
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Figure  5  2.  Filtering  of  signnl  rlt’pniricn!  noisy  images. 


by  signal-dependent  noise  modeled  by  Eq.  (5.27)  when  k  =  2.  Figure 
5-2a  is  the  original  image,  Fig.  5-2b  is  the  noisy  image  and  Fig. 
5-2c  is  the  filtered  result.  Figure  5-3a  is  a  plot  of  the  brightness 
cross-section  of  Fig.  5-2a,  and  the  subsequent  plots  are  brightness 
cross-sections  of  noisy  and  filtered  images.  The  filter  has  reduced 
the  mean  square  error  by  a  factor  of  ten. 

5. 4  Restoration  of  Noisy  Blurred  Images 

A  noisy  blurred  image  can  be  expressed  in  the  form 

£  =  H  L  +  H  (5.28) 

for  a  discrete  model,  where  f_  and  £  are  samples  of  the  object  and 
image  functions  given  by  (5.20a)  and  (5.20b),  respectively.  The 
following  fidelity  criteria  for  image  restoration  may  then  be 
formulated  as  minimization  of 

nXN  2 

j  (f"(x))  dx  (5.29) 

X1 

2 

among  all  twice  differentiable  function  feC  such  that 

I!  D_1(H  f_-£)||2  s  S  (5.30) 

where  D  is  defined  byEq.  (5.20c).  This  minimization  criteria  is 


similar  to  what  was  defined  for  unblurred  noisy  images  with  a  slight 
modification.  If  H  =  I,  the  condition  (5.  30)  reduces  to  (5.  1  3).  Using 
the  same  procedure  as  before,  the  objective  function 
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(5.  31) 


J  (f'fxD^dx  +  p[  |j_D  1  (H  f_-£)||2  +  z2  -  S] 


must  be  minimized  with  respect  to  f_.  Assuming  the  function  f  to  be 
composed  of  piecewise  polynomials  given  by  (5.  17),  the  equations 


Tc  =  Q  a 


(5.32) 


t  -2 

Qc  =  pH  D  (jr-Ha) 


(5.33) 


are  obtained,  where 


£  =  [c2,c3. 


a  =  [aj.a2, 


CN-l^ 

t 


(5.  34) 


(5. 35) 


and  Q  and  _T  are  previously  defined.  In  order  to  satisfy  condition 
(5.  18),  Cj  and  c^  are  assumed  to  be  zero.  Since  _T  is  a  positive 


definite  matrix 


c  =  r’^a 


(5.  36) 


and  substituting  (5,  36)  in  (5.  33)  yields 


Q  T'^a  =  pHt’D_2(^-H  a) 


(5.  37) 


From  (5.  17),  it  is  clear  that  a  =  f  ,  therefore 


f_  =  (HfcD’2H  +  XQ  J^QV^D"2^  (5.38) 


where  X  =  p  .  Setting  H_  =  I  leads  to  the  estimate  of  Eq.  (5.  26). 
Comparison  of  this  restoration  filter  to  the  Wiener  filter  given  by 


(5.7)  shows  the  similarity  of  these  two  filters  with  C.  replaced  by 


XQ  T  In  order  to  obtain  the  filter  numerically,  it  is  necessary 

to  invert  the  matrix  T_  as  well  as  H^D  H  +  X Q T  Q*.  Matrix  T  is  a 
positive -definite  tridiagonal  matrix  and  efficient  techniques  exist  for 
the  computation  of  its  inverse  [57]  .  Since  H  is  a  large  matrix  for 
ordinary  size  images,  and  also  H  is  not  of  full  column  rank  for  non¬ 
zero  background  pictures,  the  matrix  If  D  H  +  XR_  _T  Q*"  may  be 
ill-conditioned  or  singular  and  thus  pseudo-inversion  is  used  instead 
of  exact  inversion.  In  section  5.  6  the  experimental  results  of  this 
filter  are  presented. 

5.  5  The  Effect  of  Fidelity  Criteria  on  High  Frequency  Suppression 
In  Chapter  4,  describing  the  restoration  of  noiseless  images,  a 
fidelity  criteria  based  on  the  minimization  of  the  L^-norm  of  the 
function  was  defined  which  led  to  the  filter  given  by  (4.  19)  and 
pseudoinversion.  In  this  chapter,  the  fidelity  criterion  is  based  on 
the  minimization  of  L^-norm  of  the  second  derivative  of  the  function, 
leading  to  spline  filtering  and  restoration.  In  this  section,  the  effects 
of  different  criteria  on  the  frequency  content  of  the  restored  image 
are  discussed. 

The  objective  function  for  both  cases  has  the  one  common  term 

X  J  (f^'(x))^dx  (5.39) 

where  X  is  a  weighting  factor.  If  k  =  0,  the  derived  filter  is  given  by 
(4.  19)  or  by  pseudo-inversion,  and  k  =  2  leads  to  spline  restoration. 


r>0 


Minimization  of  the  objective  function  implies  the  reduction  of  term 
(5.  3°))  in  addition  to  the  other  terms.  Taking  the  Fourier  transform  ot 
this  term  and  applying  Parseval*s  theorem,  it  is  equivalent  to 

Xj  [(2ttu)^F(u)]  du  (5.40) 

where  u  is  the  spatial  frequency  and  F(u)  is  the  Fourier  transform  of 
fix).  As  (5.40)  shows,  the  amplitude  F(u)  at  spatial  frequency  u  is 
weighted  by  a  factor  (2ttu)  and  then  squared  and  integrated.  There¬ 
fore,  as  k  increases,  the  weighting  factor  of  high  frequency  amplitude 
increases  too.  Thus,  reducing  (5.40)  with  larger  k  implies  the 
suppression  of  the  higher  frequencies  more  than  the  lower  frequencies. 
When  k  =  0,  the  amplitudes  of  all  the  frequencies  are  equally  weighted. 

In  most  image  restoration  problems  with  white  noise,  the  original 
object  is  a  relatively  low  frequency  signal,  and  noise  has  a  flat  power 
spectral  density.  When  the  noise  power  is  zero  or  very  low,  there  is 
no  need  to  suppress  the  high  frequency  content  of  the  image.  More¬ 
over,  preserving  the  higher  frequencies  preserves  the  fine  details  of 
the  image,  and  thus,  k  =  0  is  a  good  choice  for  noiseless  images. 
Conversely,  for  noisy  images,  where  most  of  the  high  frequency 
content  of  the  image  belongs  to  the  noise,  k  =  2  is  a  good  choice. 
Generally,  as  k  increases,  the  filtered  signal  is  more  smooth  and 
correlated  with  less  higher  frequency  content. 
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5.  6  Experi  mental  Results 


In  this  section  the  experimental  results  of  restoration  of  noisy 
blurred  images  by  spline  functions  are  presented.  Figure  5 -4a 
illustrates  a  scene  of  original  size  128  y  128  pixels  enlarged  to 
256  x  256  by  cubic  spline  interpolation  for  display  purposes.  Figure 
5-4b  represents  the  image  following  motion  degradation  with  a  point- 
spread  function  of  length  7  pixels  and  the  addition  of  zero  mean 
Gaussian  noise  with  standard  deviation  1.  The  observed  image  has  a 
nonzero  background,  and  thus  the  system  is  underdetermined.  For 
restoration,  the  smoothing  parameter  X  of  the  filter  (5.  38)  is  set  to 

.01  and  .001,  and  the  SVD  algorithm  is  employed  to  compute  the 

t  -2  -It 

pseudo-inverse  of  the  matrix  H  _D  H  +  XQ  _T  Q  . 

Figure  5-5  shows  the  behavior  of  singular  values  for  different 

values  of  X  .  When  X  is  small,  it  only  affects  the  smaller  singular 

values  (right  side  portion  of  the  curves)  which  correspond  to  high 

frequency  eigenvectors,  while  for  large  X,  the  smoothing  term 

-It  t  -2 

Q  _T  Q  completely  dominates  the  deblurring  term  H  D  FI  and  its 

singular  values.  The  effect  of  X  extends  from  the  left  side  portion  of 

the  curves  to  the  right  portion  as  X  increases.  For  the  computation 

of  the  pseudoinverse  a  judicious  threshold  e  must  be  selected  for  the 

cutoff  of  singular  values.  Figures  5-4c  and  5-4d  show  the  restored 

image  for  X  =  .01,  c  =  .005,  and  for  X  =  .001,  e  =  .001,  respectively. 

As  X  and  c  decrease,  the  restored  image  becomes  sharper  with 
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(b)  Noisy  blurred,  L 


(a)  Original 


005  (d)  Restored.  X  =  .001 


Figure  5-4.  Restoration  of  noisy  blurred  images  by  spline  filter 


more  details,  but  at  the  same  time  the  amplitude  of  unwanted  high 
frequency  components  increases.  For  large  \  anc.  e,  the  restored 
image  is  smoother  and  more  correlated.  Thus,  by  choosing  these 
two  parameters  as  well  as  the  local  smoothing  window  &.,  one  can 
control  the  tradeoff  between  smoothing  and  resolution  locally  and 
globally.  It  is  obvious  that  the  quality  of  the  restored  image  is 
highly  dependent  on  the  proper  selection  of  these  parameters. 
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Chapter  6 

CONCLUSIONS  ANI)  SHOCKS  HONS  FOR 
FUR  THER  RESEARCH 

This  dissertation  has  presented  a  theoretical  and  experimental 
analysis  of  image  restoration  in  the  spline  domain.  Representing 
object,  image  and  point-spread  functions  by  splines  has  led  to  a 
more  accurate  and  realistic  model.  The  interesting  properties  of 
spline  functions,  particularly  B -splines,  have  been  used  in  image 
modeling  and  restoration. 

The  linear  integral  equation  that  describes  the  image  formation 
has  been  discretized  for  processing  by  a  digital  computer.  Various 
methods  such  as  pulse  approximation,  Newton-Cotes  and  monospline 
quadrature  formulae  are  discussed  and  compared  with  each  other. 

The  first  two  methods  are  special  and  extreme  cases  of  the  third 
one.  It  is  shown  that  the  monospline  quadrature  formulae  of  degree  3 
produces  less  error  overall  than  the  other  two  methods,  although  the 
best  choice  of  degree  for  the  monospline  is  dependent  on  the  variation 
of  the  integrand  and  sampling  interval.  For  the  class  of  rapidly 
varying  functions,  a  smaller  m  is  advised,  but  for  the  class  of 
slowly  varying  functions,  large  values  of  m  give  better  results. 

These  results  are  true  for  both  undersampled  and  oversampled 
functions, 
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B -splines  have  been  used  for  interpolation  and  approximation  of 
object  and  image  functions.  B -splines  hav*7*  some  advantages  over  the 
sine  and  cosine  basis  functions  of  the  Fourier  domain.  First,  they 
are  strictly  positive  and  thus  they  are  a  better  representation  of  image- 
intensity.  Second,  the  shape  of  these  functions  is  fixed  except  for  a 
shift  of  location.  This  property  results  in  circulant  matrices  which 
are  easy  to  handle  numerically.  Third,  their  local  basis  properties 
results  in  banded  matrices.  Fourth,  their  convolutional  property 
represents  the  convolution  integral  of  space-invariant  degradations. 

It  has  been  shown  that  using  B-splincs  and  exploiting  the  convolutional 
property  for  restoration  of  an  analytically  blurred  image  results  in  a 
more  accurate  reconstruction  of  the  original  object.  Since  the  blurred 
image  is  generally  smoother  than  the  object,  a  higher  degree  B-spline 
represents  it  better  than  the  one  representing  the  object.  The  degree 
of  B-spline  must  be  selected  according  to  the  frequency  content  (or 
variation)  of  the  approximated  function. 

The  minimum  norm  criteria  leading  to  the  pseudo-inversion  has 
been  used  for  the  restoration  of  space-variant  degradations,  over¬ 
determined  models  and  underdetermined  models.  A  numerical 
technique  called  singular-value-decomposition  is  used  for  computa¬ 
tion  of  the  pseudo-inverse.  Due  to  the  singularity  of  the  system, 
ill-conditioning  or  roundoff  error,  a  judicious  threshold  for  the 
nonzero  singular  values  must  be  made.  A  sudden  decrease  in  their 
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values  and  the  machine  precision  are  the  determining  factors.  The 
strikingly  good  results  for  the  restoration  of  astigmatism  which  is 
strongly  space -variant  show  the  validity  of  the  objective  criteria  and 
the  capability  of  the  numerical  technique. 

A  constrained  least  squares  criteria  based  on  the  minimization 
of  the  second  derivative  and  local  statistics  of  the  noise  has  been 
defined  and  the  optimal  filter,  which  produces  a  cubic  spline  function 
as  the  estimate,  has  been  derived.  This  filter  is  applicable  to  both 
space-variant  and  space-invariant  degradations.  The  parameters  of 
the  filter  determine  the  local  smoothing  window  and  overall  extent  of 
smoothing  and  thus  the  tradeoff  between  resolution  and  smoothing  is 
controllable  both  locally  and  globally.  These  properties  have  made 
the  spline  filter  capable  of  restoring  images  with  signal-dependent  or 
multiplicative  noise.  This  filter  has  been  successfully  applied  for 
filtering  images  degraded  by  film-grain  noise  which  is  modeled  as 
signal-dependent  or  multiplicative  noise.  Since  a  local  spatial 
average  has  been  used  as  a  nonstationary  estimate  of  the  signal  mean 
and  local  smoothing  factor,  the  overall  filtering  is  nonlinear.  The 
filter  has  been  also  applied  to  images  degraded  by  motion  blur  and 
noise.  The  results  show  that  as  the  smoothing  parameter  decreases, 
the  restored  image  becomes  sharper  with  more  details,  but  at  the 
same  time  the  amplitude  of  unwanted  high  frequency  components 
increases.  For  large  smoothing  parameters,  the  restored  images 
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become  more  correlated.  Finally,  the  effect  of  fidelity  criteria  on 
high  frequency  components  of  the  restored  image  has  been  analyzed. 
If  the  noise  power  is  zero  or  almost  zero,  the  fidelity  criteria  based 
on  the  minimization  of  the  function  is  advised,  but  for  noisy  images, 
the  fidelity  criteria  with  respect  to  the  minimization  of  the  second 
derivative  of  the  function  gives  better  results. 

1’he  research  pursued  in  this  dissertation  may  be  extended  in 
several  directions.  A  more  detailed  study  of  monospline  quadrature 
formulae,  particularly  free  nodes,  would  be  of  considerable  interest. 
The  number  of  nodes  and  their  locations  have  substantial  effect  on 
quadrature  error  and  accuracy  of  the  discrete  model.  This  study  can 
be  extended  with  respect  to  the  behavior  and  local  properties  of  the 
integrand.  In  approximating  the  image  function  by  splines,  variable 
sampling  may  be  studied  for  optimal  knot  placement  and  error  re¬ 
duction.  A  discussion  of  the  problem  as  well  as  some  algorithms 
may  be  found  in  references  [60]  ,  [61  ]  . 

A  more  detailed  experimental  study  of  the  spline  filter  is 
particularly  useful  in  determining  the  parameters  of  the  filter 
according  to  the  local  properties  of  the  image  and  noise  statistics. 

The  filter  can  be  applied  to  non-uniform  sampling  with  minor 
modification  of  the  filter  matrices  [Appendix  A]  .  Thus  the 
sampling  intervals  may  be  determined  with  respect  to  local  behavior 


of  the  observed  image  and  then  the  modified  filter  can  be  used. 
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Recursive  computational  techniques  are  of  special  interest  in 
real-time  processing  of  the  images.  A  version  of  the  filter  may  be 
developed  for  restoring  images  when  the  stream  of  data  comes  from 
a  scanner  or  transmitter.  Another  possibility  would  be  an  exploration 
of  the  iterative  method  for  solution  of  the  filter.  Such  a  method  has 
been  developed  for  computation  of  the  pseudoinverse  [62]  ,  [63]  . 

Since  the  matrices  of  the  spline  filter  are  banded  circulant  or 
Toeplitz,  very  efficient  numerical  techniques,  particularly  for 
space-invariant  degradations,  may  be  used. 

As  mentioned  earlier,  spline  functions  have  many  desirable 
approximating  and  interpolating  characteristics  that  can  be  used  in 
various  areas  of  digital  image  processing.  So  far,  some  attention 
has  been  given  to  potential  applications  of  spline  functions  in  image 
restoration  and  enlargement  [45]  ,  [64]  .  Although  this  is  just 
beginning  work  in  utilizing  the  properties  of  spline  functions  in  image 
restoration,  the  successful  results  indicate  that  further  fruitful 
research  can  be  performed  in  this  field.  Spline  functions  may  be  used 
in  other  areas  of  digital  image  processing  such  as  image  coding  and 


reconstruction  of  images  from  their  projections. 


APPENDIX  A 


1 


Using  vivtnr  space  notation,  Fqs.  (A-7)  and  (A-8)  can  be  written  as 

_Tc  =  o'a  IA-9) 

-2 

Q  c_  =  pIO  (g_-a)  (A-10) 


where  £  =  [c^.  .  .  .  ,  c  a  =  [a  ^ . a^]*",  _T_  is  a  positive  de¬ 

finite  tridiagonal  matrix  of  order  N-2  with  elements 


t  . 
1 1 


and  Q  is  a  N  X  (N-2)  tridiagonal  matrix  with  elements 

1  1  _j_  1 

qi,i  +  l  =  h.  ’  qi+l,i+l  "  h.  .  "  h.  ’  V2,i  +  1  "  h.  ,  * 

i  i +l  i  i+t 

t  2 

Premultiplying  both  sides  of  Eq.  (A-10)  by  Q  D 

Q 1 10  Q_£  =  pQ^  -  pQ*a  ( A  - 1 1 ) 

t 

and  substituting  CO  a  from  (A  -  9)  into  (A  - 1 1 )  gives 

(0*13  0  +  p_T)£  =  pQ l£  •  (A- 12) 

t  2 

Since  Q  10  Q  +  p_T  is  positive  definite,  c_  amy  be  obtained  as 

£  =  p(QtlD^Q  +  pT)  *Q*£  .  (A-14) 

From  (A-10),  a_  is  derived  as 

-1  2 

a  =  £  -  p  D  Q  c  (A-14) 

Substituting  ( A  -  1 3 )  in  (A-14)  gives 
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APPENDIX  B 


MATRIX  IDENTITY 

Lemma:  For  any  nonsingular  matrices  _B  and  £,  the  following 
identity 

-1  -1  t  -1  t  -It 

(C  +A  B_  A  )  =  £  -  £  A  (  A  £  A +B_I  A  £ 

-1  -1  t  -1  t  -1 

is  valid  for  any  matrix  A  if  (£  +  A  B_  A  )  and  ( A^  £  A+B_) 

Proof:  Assuming  an  auxiliary  matrix  13  of  the  form 

-1  -It 

D  =  £  +  A  B  A 

and  premultiplying  both  sides  of  Eq.  (B-l)  by  D  gives 

-1  -1  -1  -1  t 

1  =  D  C  +  D  A_B  A  . 

Postmultiplying  (B-2)  by  C^,  the  equation 

-1  -1  -1  t 

C  =  D  +  D  A  B_  A  £ 

is  obtained,  and  this  is  rewritten  as 

D_1  =  £  -  D ' 1 A  1 A lC  . 

Postmultiplying  (B-3)  by  A_  gives 

-1  -1  -1  t 

£  A  =  D  A  +  D  A  £  A  £  A 

-1  -It 

=  D  A  (I  +  B  A  £  A  ) 

=  R  1 A  B  \b +Al£  A)  . 


matrix 


exist. 


(B-l) 


(B-2 1 


(B-3) 


(B-4) 


(B-5) 
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Since  the  inverse  ofB_+A^£Aexists,  D_  £  B_  is  obtai  ned  a  s 

_D~ 1 A  B  =  £  A(B_  +  AC_  A)’1  (B-6) 

and  this  can  be  substituted  in  (B-4)  to  give 

D"1  =  £_-  £  A(B_  +  VcAf'/c  (B-7) 

or  from  (B  - 1 ) 

(C_1  +  A  =  £  -  £  A(B  +At£A)'1AtC  .  (B-8) 

Q.  E.D. 
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